home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / integration / test.c < prev    next >
Encoding:
C/C++ Source or Header  |  2001-08-22  |  75.8 KB  |  2,106 lines
  1. /* integration/test.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21. #include <stdlib.h>
  22. #include <stdio.h>
  23. #include <math.h>
  24. #include <gsl/gsl_math.h>
  25. #include <gsl/gsl_integration.h>
  26. #include <gsl/gsl_errno.h>
  27. #include <gsl/gsl_test.h>
  28. #include <gsl/gsl_ieee_utils.h>
  29.  
  30. #include "tests.h"
  31.  
  32. gsl_function make_function (double (* f) (double, void *), double * p);
  33.  
  34. gsl_function make_function (double (* f) (double, void *), double * p)
  35. {
  36.   gsl_function f_new;
  37.  
  38.   f_new.function = f ;
  39.   f_new.params = p ;
  40.  
  41.   return f_new;
  42. }
  43.  
  44. struct counter_params {
  45.   gsl_function * f;
  46.   int neval;
  47. } ;
  48.  
  49. double counter (double x, void * params);
  50. gsl_function make_counter (gsl_function * f, struct counter_params * p);
  51.  
  52. double 
  53. counter (double x, void * params)
  54. {
  55.   struct counter_params * p = (struct counter_params *) params;
  56.   p->neval++ ; /* increment counter */
  57.   return GSL_FN_EVAL(p->f, x);
  58. }
  59.  
  60. gsl_function make_counter (gsl_function * f, struct counter_params * p)
  61. {
  62.   gsl_function f_new;
  63.  
  64.   p->f = f;
  65.   p->neval = 0 ;
  66.   
  67.   f_new.function = &counter ;
  68.   f_new.params = p ;
  69.  
  70.   return f_new;
  71. }
  72.  
  73. void my_error_handler (const char *reason, const char *file,
  74.                int line, int err);
  75.  
  76. int main (void)
  77. {
  78.   gsl_ieee_env_setup ();
  79.   gsl_set_error_handler (&my_error_handler); 
  80.  
  81.   /* Test the basic Gauss-Kronrod rules with a smooth positive function. */
  82.  
  83.   {
  84.     double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
  85.     double exp_result = 7.716049357767090777E-02;
  86.     double exp_abserr = 2.990224871000550874E-06;
  87.     double exp_resabs = 7.716049357767090777E-02;
  88.     double exp_resasc = 4.434273814139995384E-02;
  89.  
  90.     double alpha = 2.6 ;
  91.     gsl_function f = make_function(&f1, &alpha) ;
  92.  
  93.     gsl_integration_qk15 (&f, 0.0, 1.0, 
  94.                   &result, &abserr, &resabs, &resasc) ;
  95.     gsl_test_rel(result,exp_result,1e-15,"qk15(f1) smooth result") ;
  96.     gsl_test_rel(abserr,exp_abserr,1e-7,"qk15(f1) smooth abserr") ;
  97.     gsl_test_rel(resabs,exp_resabs,1e-15,"qk15(f1) smooth resabs") ;    
  98.     gsl_test_rel(resasc,exp_resasc,1e-15,"qk15(f1) smooth resasc") ;
  99.  
  100.     gsl_integration_qk15 (&f, 1.0, 0.0, 
  101.                   &result, &abserr, &resabs, &resasc) ;
  102.  
  103.     gsl_test_rel(result,-exp_result,1e-15,"qk15(f1) reverse result") ;
  104.     gsl_test_rel(abserr,exp_abserr,1e-7,"qk15(f1) reverse abserr") ;
  105.     gsl_test_rel(resabs,exp_resabs,1e-15,"qk15(f1) reverse resabs") ;    
  106.     gsl_test_rel(resasc,exp_resasc,1e-15,"qk15(f1) reverse resasc") ;
  107.   }
  108.  
  109.   {
  110.     double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
  111.     double exp_result = 7.716049379303084599E-02;
  112.     double exp_abserr = 9.424302194248481445E-08;
  113.     double exp_resabs = 7.716049379303084599E-02;
  114.     double exp_resasc = 4.434311425038358484E-02;
  115.  
  116.     double alpha = 2.6 ;
  117.     gsl_function f = make_function(&f1, &alpha);
  118.  
  119.     gsl_integration_qk21 (&f, 0.0, 1.0, 
  120.                   &result, &abserr, &resabs, &resasc) ;
  121.     gsl_test_rel(result,exp_result,1e-15,"qk21(f1) smooth result") ;
  122.     gsl_test_rel(abserr,exp_abserr,1e-7,"qk21(f1) smooth abserr") ;
  123.     gsl_test_rel(resabs,exp_resabs,1e-15,"qk21(f1) smooth resabs") ;    
  124.     gsl_test_rel(resasc,exp_resasc,1e-15,"qk21(f1) smooth resasc") ;
  125.  
  126.     gsl_integration_qk21 (&f, 1.0, 0.0, 
  127.                   &result, &abserr, &resabs, &resasc) ;
  128.     gsl_test_rel(result,-exp_result,1e-15,"qk21(f1) reverse result") ;
  129.     gsl_test_rel(abserr,exp_abserr,1e-7,"qk21(f1) reverse abserr") ;
  130.     gsl_test_rel(resabs,exp_resabs,1e-15,"qk21(f1) reverse resabs") ;    
  131.     gsl_test_rel(resasc,exp_resasc,1e-15,"qk21(f1) reverse resasc") ;
  132.   }
  133.  
  134.   {
  135.     double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
  136.     double exp_result = 7.716049382494900855E-02;
  137.     double exp_abserr = 1.713503193600029893E-09;
  138.     double exp_resabs = 7.716049382494900855E-02;
  139.     double exp_resasc = 4.427995051868838933E-02;
  140.  
  141.     double alpha = 2.6 ;
  142.     gsl_function f = make_function(&f1, &alpha);
  143.  
  144.     gsl_integration_qk31 (&f, 0.0, 1.0, 
  145.                   &result, &abserr, &resabs, &resasc) ;
  146.     gsl_test_rel(result,exp_result,1e-15,"qk31(f1) smooth result") ;
  147.     gsl_test_rel(abserr,exp_abserr,1e-7,"qk31(f1) smooth abserr") ;
  148.     gsl_test_rel(resabs,exp_resabs,1e-15,"qk31(f1) smooth resabs") ;    
  149.     gsl_test_rel(resasc,exp_resasc,1e-15,"qk31(f1) smooth resasc") ;
  150.  
  151.     gsl_integration_qk31 (&f, 1.0, 0.0, 
  152.                   &result, &abserr, &resabs, &resasc) ;
  153.     gsl_test_rel(result,-exp_result,1e-15,"qk31(f1) reverse result") ;
  154.     gsl_test_rel(abserr,exp_abserr,1e-7,"qk31(f1) reverse abserr") ;
  155.     gsl_test_rel(resabs,exp_resabs,1e-15,"qk31(f1) reverse resabs") ;    
  156.     gsl_test_rel(resasc,exp_resasc,1e-15,"qk31(f1) reverse resasc") ;
  157.   }
  158.  
  159.   {
  160.     double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
  161.     double exp_result = 7.716049382681375302E-02;
  162.     double exp_abserr = 9.576386660975511224E-11;
  163.     double exp_resabs = 7.716049382681375302E-02;
  164.     double exp_resasc = 4.421521169637691873E-02;
  165.  
  166.     double alpha = 2.6 ;
  167.     gsl_function f = make_function(&f1, &alpha);
  168.  
  169.     gsl_integration_qk41 (&f, 0.0, 1.0, 
  170.                   &result, &abserr, &resabs, &resasc) ;
  171.     gsl_test_rel(result,exp_result,1e-15,"qk41(f1) smooth result") ;
  172.     gsl_test_rel(abserr,exp_abserr,1e-7,"qk41(f1) smooth abserr") ;
  173.     gsl_test_rel(resabs,exp_resabs,1e-15,"qk41(f1) smooth resabs") ;    
  174.     gsl_test_rel(resasc,exp_resasc,1e-15,"qk41(f1) smooth resasc") ;
  175.  
  176.     gsl_integration_qk41 (&f, 1.0, 0.0, 
  177.                   &result, &abserr, &resabs, &resasc) ;
  178.     gsl_test_rel(result,-exp_result,1e-15,"qk41(f1) reverse result") ;
  179.     gsl_test_rel(abserr,exp_abserr,1e-7,"qk41(f1) reverse abserr") ;
  180.     gsl_test_rel(resabs,exp_resabs,1e-15,"qk41(f1) reverse resabs") ;    
  181.     gsl_test_rel(resasc,exp_resasc,1e-15,"qk41(f1) reverse resasc") ;
  182.   }
  183.  
  184.   {
  185.     double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
  186.     double exp_result = 7.716049382708510540E-02;
  187.     double exp_abserr = 1.002079980317363772E-11;
  188.     double exp_resabs = 7.716049382708510540E-02;
  189.     double exp_resasc = 4.416474291216854892E-02;
  190.  
  191.     double alpha = 2.6 ;
  192.     gsl_function f = make_function(&f1, &alpha);
  193.  
  194.     gsl_integration_qk51 (&f, 0.0, 1.0, 
  195.                   &result, &abserr, &resabs, &resasc) ;
  196.     gsl_test_rel(result,exp_result,1e-15,"qk51(f1) smooth result") ;
  197.     gsl_test_rel(abserr,exp_abserr,1e-5,"qk51(f1) smooth abserr") ;
  198.     gsl_test_rel(resabs,exp_resabs,1e-15,"qk51(f1) smooth resabs") ;    
  199.     gsl_test_rel(resasc,exp_resasc,1e-15,"qk51(f1) smooth resasc") ;
  200.  
  201.     gsl_integration_qk51 (&f, 1.0, 0.0, 
  202.                   &result, &abserr, &resabs, &resasc) ;
  203.     gsl_test_rel(result,-exp_result,1e-15,"qk51(f1) reverse result") ;
  204.     gsl_test_rel(abserr,exp_abserr,1e-5,"qk51(f1) reverse abserr") ;
  205.     gsl_test_rel(resabs,exp_resabs,1e-15,"qk51(f1) reverse resabs") ;    
  206.     gsl_test_rel(resasc,exp_resasc,1e-15,"qk51(f1) reverse resasc") ;
  207.   }
  208.  
  209.   {
  210.     double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
  211.     double exp_result = 7.716049382713800753E-02;
  212.     double exp_abserr = 1.566060362296155616E-12;
  213.     double exp_resabs = 7.716049382713800753E-02;
  214.     double exp_resasc = 4.419287685934316506E-02;
  215.  
  216.     double alpha = 2.6 ;
  217.     gsl_function f = make_function(&f1, &alpha);
  218.  
  219.     gsl_integration_qk61 (&f, 0.0, 1.0, 
  220.                   &result, &abserr, &resabs, &resasc) ;
  221.     gsl_test_rel(result,exp_result,1e-15,"qk61(f1) smooth result") ;
  222.     gsl_test_rel(abserr,exp_abserr,1e-5,"qk61(f1) smooth abserr") ;
  223.     gsl_test_rel(resabs,exp_resabs,1e-15,"qk61(f1) smooth resabs") ;    
  224.     gsl_test_rel(resasc,exp_resasc,1e-15,"qk61(f1) smooth resasc") ;
  225.  
  226.     gsl_integration_qk61 (&f, 1.0, 0.0, 
  227.                   &result, &abserr, &resabs, &resasc) ;
  228.     gsl_test_rel(result,-exp_result,1e-15,"qk61(f1) reverse result") ;
  229.     gsl_test_rel(abserr,exp_abserr,1e-5,"qk61(f1) reverse abserr") ;
  230.     gsl_test_rel(resabs,exp_resabs,1e-15,"qk61(f1) reverse resabs") ;    
  231.     gsl_test_rel(resasc,exp_resasc,1e-15,"qk61(f1) reverse resasc") ;
  232.   }
  233.  
  234.   /* Now test the basic rules with a positive function that has a
  235.      singularity. This should give large values of abserr which would
  236.      find discrepancies in the abserr calculation. */
  237.  
  238.   {
  239.     double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
  240.     double exp_result = 1.555688196612745777E+01;
  241.     double exp_abserr = 2.350164577239293706E+01;
  242.     double exp_resabs = 1.555688196612745777E+01;
  243.     double exp_resasc = 2.350164577239293706E+01;
  244.  
  245.     double alpha = -0.9 ;
  246.     gsl_function f = make_function(&f1, &alpha);
  247.  
  248.     gsl_integration_qk15 (&f, 0.0, 1.0, 
  249.                   &result, &abserr, &resabs, &resasc) ;
  250.     gsl_test_rel(result,exp_result,1e-15,"qk15(f1) singular result") ;
  251.     gsl_test_rel(abserr,exp_abserr,1e-7,"qk15(f1) singular abserr") ;
  252.     gsl_test_rel(resabs,exp_resabs,1e-15,"qk15(f1) singular resabs") ;    
  253.     gsl_test_rel(resasc,exp_resasc,1e-15,"qk15(f1) singular resasc") ;
  254.  
  255.     gsl_integration_qk15 (&f, 1.0, 0.0, 
  256.                   &result, &abserr, &resabs, &resasc) ;
  257.     gsl_test_rel(result,-exp_result,1e-15,"qk15(f1) reverse result") ;
  258.     gsl_test_rel(abserr,exp_abserr,1e-7,"qk15(f1) reverse abserr") ;
  259.     gsl_test_rel(resabs,exp_resabs,1e-15,"qk15(f1) reverse resabs") ;    
  260.     gsl_test_rel(resasc,exp_resasc,1e-15,"qk15(f1) reverse resasc") ;
  261.   }
  262.  
  263.   {
  264.     double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
  265.     double exp_result = 1.799045317938126232E+01;
  266.     double exp_abserr = 2.782360287710622515E+01;
  267.     double exp_resabs = 1.799045317938126232E+01;
  268.     double exp_resasc = 2.782360287710622515E+01;
  269.  
  270.     double alpha = -0.9 ;
  271.     gsl_function f = make_function(&f1, &alpha);
  272.  
  273.     gsl_integration_qk21 (&f, 0.0, 1.0, 
  274.                   &result, &abserr, &resabs, &resasc) ;
  275.     gsl_test_rel(result,exp_result,1e-15,"qk21(f1) singular result") ;
  276.     gsl_test_rel(abserr,exp_abserr,1e-7,"qk21(f1) singular abserr") ;
  277.     gsl_test_rel(resabs,exp_resabs,1e-15,"qk21(f1) singular resabs") ;    
  278.     gsl_test_rel(resasc,exp_resasc,1e-15,"qk21(f1) singular resasc") ;
  279.  
  280.     gsl_integration_qk21 (&f, 1.0, 0.0, 
  281.                   &result, &abserr, &resabs, &resasc) ;
  282.     gsl_test_rel(result,-exp_result,1e-15,"qk21(f1) reverse result") ;
  283.     gsl_test_rel(abserr,exp_abserr,1e-7,"qk21(f1) reverse abserr") ;
  284.     gsl_test_rel(resabs,exp_resabs,1e-15,"qk21(f1) reverse resabs") ;    
  285.     gsl_test_rel(resasc,exp_resasc,1e-15,"qk21(f1) reverse resasc") ;
  286.   }
  287.  
  288.   {
  289.     double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
  290.     double exp_result = 2.081873305159121657E+01;
  291.     double exp_abserr = 3.296500137482590276E+01;
  292.     double exp_resabs = 2.081873305159121301E+01;
  293.     double exp_resasc = 3.296500137482590276E+01;
  294.  
  295.     double alpha = -0.9 ;
  296.     gsl_function f = make_function(&f1, &alpha);
  297.  
  298.     gsl_integration_qk31 (&f, 0.0, 1.0, 
  299.                   &result, &abserr, &resabs, &resasc) ;
  300.     gsl_test_rel(result,exp_result,1e-15,"qk31(f1) singular result") ;
  301.     gsl_test_rel(abserr,exp_abserr,1e-7,"qk31(f1) singular abserr") ;
  302.     gsl_test_rel(resabs,exp_resabs,1e-15,"qk31(f1) singular resabs") ;    
  303.     gsl_test_rel(resasc,exp_resasc,1e-15,"qk31(f1) singular resasc") ;
  304.  
  305.     gsl_integration_qk31 (&f, 1.0, 0.0, 
  306.                   &result, &abserr, &resabs, &resasc) ;
  307.     gsl_test_rel(result,-exp_result,1e-15,"qk31(f1) reverse result") ;
  308.     gsl_test_rel(abserr,exp_abserr,1e-7,"qk31(f1) reverse abserr") ;
  309.     gsl_test_rel(resabs,exp_resabs,1e-15,"qk31(f1) reverse resabs") ;    
  310.     gsl_test_rel(resasc,exp_resasc,1e-15,"qk31(f1) reverse resasc") ;
  311.   }
  312.  
  313.   {
  314.     double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
  315.     double exp_result = 2.288677623903126701E+01;
  316.     double exp_abserr = 3.671538820274916048E+01;
  317.     double exp_resabs = 2.288677623903126701E+01;
  318.     double exp_resasc = 3.671538820274916048E+01;
  319.  
  320.     double alpha = -0.9 ;
  321.     gsl_function f = make_function(&f1, &alpha);
  322.  
  323.     gsl_integration_qk41 (&f, 0.0, 1.0, 
  324.                   &result, &abserr, &resabs, &resasc) ;
  325.     gsl_test_rel(result,exp_result,1e-15,"qk41(f1) singular result") ;
  326.     gsl_test_rel(abserr,exp_abserr,1e-7,"qk41(f1) singular abserr") ;
  327.     gsl_test_rel(resabs,exp_resabs,1e-15,"qk41(f1) singular resabs") ;    
  328.     gsl_test_rel(resasc,exp_resasc,1e-15,"qk41(f1) singular resasc") ;
  329.  
  330.     gsl_integration_qk41 (&f, 1.0, 0.0, 
  331.                   &result, &abserr, &resabs, &resasc) ;
  332.     gsl_test_rel(result,-exp_result,1e-15,"qk41(f1) reverse result") ;
  333.     gsl_test_rel(abserr,exp_abserr,1e-7,"qk41(f1) reverse abserr") ;
  334.     gsl_test_rel(resabs,exp_resabs,1e-15,"qk41(f1) reverse resabs") ;    
  335.     gsl_test_rel(resasc,exp_resasc,1e-15,"qk41(f1) reverse resasc") ;
  336.   }
  337.  
  338.   {
  339.     double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
  340.     double exp_result = 2.449953612016972215E+01;
  341.     double exp_abserr = 3.967771249391228849E+01;
  342.     double exp_resabs = 2.449953612016972215E+01;
  343.     double exp_resasc = 3.967771249391228849E+01;
  344.  
  345.     double alpha = -0.9 ;
  346.     gsl_function f = make_function(&f1, &alpha);
  347.  
  348.     gsl_integration_qk51 (&f, 0.0, 1.0, 
  349.                   &result, &abserr, &resabs, &resasc) ;
  350.     gsl_test_rel(result,exp_result,1e-15,"qk51(f1) singular result") ;
  351.     gsl_test_rel(abserr,exp_abserr,1e-7,"qk51(f1) singular abserr") ;
  352.     gsl_test_rel(resabs,exp_resabs,1e-15,"qk51(f1) singular resabs") ;    
  353.     gsl_test_rel(resasc,exp_resasc,1e-15,"qk51(f1) singular resasc") ;
  354.  
  355.     gsl_integration_qk51 (&f, 1.0, 0.0, 
  356.                   &result, &abserr, &resabs, &resasc) ;
  357.     gsl_test_rel(result,-exp_result,1e-15,"qk51(f1) reverse result") ;
  358.     gsl_test_rel(abserr,exp_abserr,1e-7,"qk51(f1) reverse abserr") ;
  359.     gsl_test_rel(resabs,exp_resabs,1e-15,"qk51(f1) reverse resabs") ;    
  360.     gsl_test_rel(resasc,exp_resasc,1e-15,"qk51(f1) reverse resasc") ;
  361.   }
  362.  
  363.   {
  364.     double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
  365.     double exp_result = 2.583030240976628988E+01;
  366.     double exp_abserr = 4.213750493076978643E+01;
  367.     double exp_resabs = 2.583030240976628988E+01;
  368.     double exp_resasc = 4.213750493076978643E+01;
  369.  
  370.     double alpha = -0.9 ;
  371.     gsl_function f = make_function(&f1, &alpha);
  372.  
  373.     gsl_integration_qk61 (&f, 0.0, 1.0, 
  374.                   &result, &abserr, &resabs, &resasc) ;
  375.     gsl_test_rel(result,exp_result,1e-15,"qk61(f1) singular result") ;
  376.     gsl_test_rel(abserr,exp_abserr,1e-7,"qk61(f1) singular abserr") ;
  377.     gsl_test_rel(resabs,exp_resabs,1e-15,"qk61(f1) singular resabs") ;    
  378.     gsl_test_rel(resasc,exp_resasc,1e-15,"qk61(f1) singular resasc") ;
  379.  
  380.     gsl_integration_qk61 (&f, 1.0, 0.0, 
  381.                   &result, &abserr, &resabs, &resasc) ;
  382.     gsl_test_rel(result,-exp_result,1e-15,"qk61(f1) reverse result") ;
  383.     gsl_test_rel(abserr,exp_abserr,1e-7,"qk61(f1) reverse abserr") ;
  384.     gsl_test_rel(resabs,exp_resabs,1e-15,"qk61(f1) reverse resabs") ;    
  385.     gsl_test_rel(resasc,exp_resasc,1e-15,"qk61(f1) reverse resasc") ;
  386.   }
  387.  
  388.   /* Test the basic Gauss-Kronrod rules with a smooth oscillating
  389.      function, over an unsymmetric range. This should find any
  390.      discrepancies in the abscissae. */
  391.  
  392.   {
  393.     double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
  394.     double exp_result =-7.238969575483799046E-01;
  395.     double exp_abserr = 8.760080200939757174E-06;
  396.     double exp_resabs = 1.165564172429140788E+00;
  397.     double exp_resasc = 9.334560307787327371E-01;
  398.  
  399.     double alpha = 1.3 ;
  400.     gsl_function f = make_function(&f3, &alpha);
  401.  
  402.     gsl_integration_qk15 (&f, 0.3, 2.71, 
  403.                   &result, &abserr, &resabs, &resasc) ;
  404.     gsl_test_rel(result,exp_result,1e-15,"qk15(f3) oscill result") ;
  405.     gsl_test_rel(abserr,exp_abserr,1e-7,"qk15(f3) oscill abserr") ;
  406.     gsl_test_rel(resabs,exp_resabs,1e-15,"qk15(f3) oscill resabs") ;
  407.     gsl_test_rel(resasc,exp_resasc,1e-15,"qk15(f3) oscill resasc") ;
  408.  
  409.     gsl_integration_qk15 (&f, 2.71, 0.3, 
  410.                   &result, &abserr, &resabs, &resasc) ;
  411.     gsl_test_rel(result,-exp_result,1e-15,"qk15(f3) reverse result") ;
  412.     gsl_test_rel(abserr,exp_abserr,1e-7,"qk15(f3) reverse abserr") ;
  413.     gsl_test_rel(resabs,exp_resabs,1e-15,"qk15(f3) reverse resabs") ;
  414.     gsl_test_rel(resasc,exp_resasc,1e-15,"qk15(f3) reverse resasc") ;
  415.   }
  416.  
  417.   {
  418.     double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
  419.     double exp_result =-7.238969575482959717E-01;
  420.     double exp_abserr = 7.999213141433641888E-11;
  421.     double exp_resabs = 1.150829032708484023E+00;
  422.     double exp_resasc = 9.297591249133687619E-01;
  423.  
  424.     double alpha = 1.3 ;
  425.     gsl_function f = make_function(&f3, &alpha);
  426.     
  427.     gsl_integration_qk21 (&f, 0.3, 2.71, 
  428.                   &result, &abserr, &resabs, &resasc) ;
  429.     gsl_test_rel(result,exp_result,1e-15,"qk21(f3) oscill result") ;
  430.     gsl_test_rel(abserr,exp_abserr,1e-5,"qk21(f3) oscill abserr") ;
  431.     gsl_test_rel(resabs,exp_resabs,1e-15,"qk21(f3) oscill resabs") ;
  432.     gsl_test_rel(resasc,exp_resasc,1e-15,"qk21(f3) oscill resasc") ;
  433.  
  434.     gsl_integration_qk21 (&f, 2.71, 0.3,
  435.                   &result, &abserr, &resabs, &resasc) ;
  436.     gsl_test_rel(result,-exp_result,1e-15,"qk21(f3) reverse result") ;
  437.     gsl_test_rel(abserr,exp_abserr,1e-5,"qk21(f3) reverse abserr") ;
  438.     gsl_test_rel(resabs,exp_resabs,1e-15,"qk21(f3) reverse resabs") ;
  439.     gsl_test_rel(resasc,exp_resasc,1e-15,"qk21(f3) reverse resasc") ;
  440.   }
  441.  
  442.   {
  443.     double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
  444.     double exp_result =-7.238969575482959717E-01;
  445.     double exp_abserr = 1.285805464427459261E-14;
  446.     double exp_resabs = 1.158150602093290571E+00;
  447.     double exp_resasc = 9.277828092501518853E-01;
  448.  
  449.     double alpha = 1.3 ;
  450.     gsl_function f = make_function(&f3, &alpha);
  451.  
  452.     gsl_integration_qk31 (&f, 0.3, 2.71, 
  453.                   &result, &abserr, &resabs, &resasc) ;
  454.     gsl_test_rel(result,exp_result,1e-15,"qk31(f3) oscill result") ;
  455.     gsl_test_rel(abserr,exp_abserr,1e-7,"qk31(f3) oscill abserr") ;
  456.     gsl_test_rel(resabs,exp_resabs,1e-15,"qk31(f3) oscill resabs") ;
  457.     gsl_test_rel(resasc,exp_resasc,1e-15,"qk31(f3) oscill resasc") ;
  458.  
  459.     gsl_integration_qk31 (&f, 2.71, 0.3, 
  460.                   &result, &abserr, &resabs, &resasc) ;
  461.     gsl_test_rel(result,-exp_result,1e-15,"qk31(f3) reverse result") ;
  462.     gsl_test_rel(abserr,exp_abserr,1e-7,"qk31(f3) reverse abserr") ;
  463.     gsl_test_rel(resabs,exp_resabs,1e-15,"qk31(f3) reverse resabs") ;
  464.     gsl_test_rel(resasc,exp_resasc,1e-15,"qk31(f3) reverse resasc") ;
  465.   }
  466.  
  467.   {
  468.     double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
  469.     double exp_result =-7.238969575482959717E-01;
  470.     double exp_abserr = 1.286535726271015626E-14;
  471.     double exp_resabs = 1.158808363486595328E+00;
  472.     double exp_resasc = 9.264382258645686985E-01;
  473.  
  474.     double alpha = 1.3 ;
  475.     gsl_function f = make_function(&f3, &alpha);
  476.  
  477.     gsl_integration_qk41 (&f, 0.3, 2.71, 
  478.                   &result, &abserr, &resabs, &resasc) ;
  479.     gsl_test_rel(result,exp_result,1e-15,"qk41(f3) oscill result") ;
  480.     gsl_test_rel(abserr,exp_abserr,1e-7,"qk41(f3) oscill abserr") ;
  481.     gsl_test_rel(resabs,exp_resabs,1e-15,"qk41(f3) oscill resabs") ;
  482.     gsl_test_rel(resasc,exp_resasc,1e-15,"qk41(f3) oscill resasc") ;
  483.  
  484.     gsl_integration_qk41 (&f, 2.71, 0.3,
  485.                   &result, &abserr, &resabs, &resasc) ;
  486.     gsl_test_rel(result,-exp_result,1e-15,"qk41(f3) reverse result") ;
  487.     gsl_test_rel(abserr,exp_abserr,1e-7,"qk41(f3) reverse abserr") ;
  488.     gsl_test_rel(resabs,exp_resabs,1e-15,"qk41(f3) reverse resabs") ;
  489.     gsl_test_rel(resasc,exp_resasc,1e-15,"qk41(f3) reverse resasc") ;
  490.   }
  491.  
  492.   {
  493.     double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
  494.     double exp_result =-7.238969575482961938E-01;
  495.     double exp_abserr = 1.285290995039385778E-14;
  496.     double exp_resabs = 1.157687209264406381E+00;
  497.     double exp_resasc = 9.264666884071264263E-01;
  498.  
  499.     double alpha = 1.3 ;
  500.     gsl_function f = make_function(&f3, &alpha);
  501.  
  502.     gsl_integration_qk51 (&f, 0.3, 2.71, 
  503.                   &result, &abserr, &resabs, &resasc) ;
  504.     gsl_test_rel(result,exp_result,1e-15,"qk51(f3) oscill result") ;
  505.     gsl_test_rel(abserr,exp_abserr,1e-7,"qk51(f3) oscill abserr") ;
  506.     gsl_test_rel(resabs,exp_resabs,1e-15,"qk51(f3) oscill resabs") ;
  507.     gsl_test_rel(resasc,exp_resasc,1e-15,"qk51(f3) oscill resasc") ;
  508.  
  509.     gsl_integration_qk51 (&f, 2.71, 0.3,
  510.                   &result, &abserr, &resabs, &resasc) ;
  511.     gsl_test_rel(result,-exp_result,1e-15,"qk51(f3) reverse result") ;
  512.     gsl_test_rel(abserr,exp_abserr,1e-7,"qk51(f3) reverse abserr") ;
  513.     gsl_test_rel(resabs,exp_resabs,1e-15,"qk51(f3) reverse resabs") ;
  514.     gsl_test_rel(resasc,exp_resasc,1e-15,"qk51(f3) reverse resasc") ;
  515.   }
  516.  
  517.   {
  518.     double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
  519.     double exp_result =-7.238969575482959717E-01;
  520.     double exp_abserr = 1.286438572027470736E-14;
  521.     double exp_resabs = 1.158720854723590099E+00;
  522.     double exp_resasc = 9.270469641771273972E-01;
  523.  
  524.     double alpha = 1.3 ;
  525.     gsl_function f = make_function(&f3, &alpha);
  526.  
  527.     gsl_integration_qk61 (&f, 0.3, 2.71, 
  528.                   &result, &abserr, &resabs, &resasc) ;
  529.     gsl_test_rel(result,exp_result,1e-15,"qk61(f3) oscill result") ;
  530.     gsl_test_rel(abserr,exp_abserr,1e-7,"qk61(f3) oscill abserr") ;
  531.     gsl_test_rel(resabs,exp_resabs,1e-15,"qk61(f3) oscill resabs") ;
  532.     gsl_test_rel(resasc,exp_resasc,1e-15,"qk61(f3) oscill resasc") ;
  533.  
  534.     gsl_integration_qk61 (&f, 2.71, 0.3,
  535.                   &result, &abserr, &resabs, &resasc) ;
  536.     gsl_test_rel(result,-exp_result,1e-15,"qk61(f3) reverse result") ;
  537.     gsl_test_rel(abserr,exp_abserr,1e-7,"qk61(f3) reverse abserr") ;
  538.     gsl_test_rel(resabs,exp_resabs,1e-15,"qk61(f3) reverse resabs") ;
  539.     gsl_test_rel(resasc,exp_resasc,1e-15,"qk61(f3) reverse resasc") ;
  540.   }
  541.  
  542.   /* Test the non-adaptive gaussian integrator QNG */
  543.  
  544.   {
  545.     int status = 0; size_t neval = 0 ;
  546.     double result = 0, abserr = 0 ;
  547.     double exp_result = 7.716049379303083211E-02;
  548.     double exp_abserr = 9.424302199601294244E-08;
  549.     int exp_neval  =  21;
  550.     int exp_ier    =   0;
  551.  
  552.     double alpha = 2.6 ;
  553.     gsl_function f = make_function(&f1, &alpha);
  554.     
  555.     status = gsl_integration_qng (&f, 0.0, 1.0, 1e-1, 0.0,
  556.                   &result, &abserr, &neval) ;
  557.     gsl_test_rel(result,exp_result,1e-15,"qng(f1) smooth result") ;
  558.     gsl_test_rel(abserr,exp_abserr,1e-7,"qng(f1) smooth abserr") ;
  559.     gsl_test_int((int)neval,exp_neval,"qng(f1) smooth neval") ;  
  560.     gsl_test_int(status,exp_ier,"qng(f1) smooth status") ;
  561.  
  562.     status = gsl_integration_qng (&f, 1.0, 0.0, 1e-1, 0.0,
  563.                   &result, &abserr, &neval) ;
  564.     gsl_test_rel(result,-exp_result,1e-15,"qng(f1) reverse result") ;
  565.     gsl_test_rel(abserr,exp_abserr,1e-7,"qng(f1) reverse abserr") ;
  566.     gsl_test_int((int)neval,exp_neval,"qng(f1) reverse neval") ;  
  567.     gsl_test_int(status,exp_ier,"qng(f1) reverse status") ;
  568.   }
  569.  
  570.   {
  571.     int status = 0; size_t neval = 0 ;
  572.     double result = 0, abserr = 0 ;
  573.  
  574.     double exp_result = 7.716049382706505200E-02;
  575.     double exp_abserr = 2.666893044866214501E-12;
  576.     int exp_neval  =  43;
  577.     int exp_ier    =   0;
  578.  
  579.     double alpha = 2.6 ;
  580.     gsl_function f = make_function(&f1, &alpha);
  581.  
  582.     status = gsl_integration_qng (&f, 0.0, 1.0, 0.0, 1e-9,
  583.                   &result, &abserr, &neval) ;
  584.     gsl_test_rel(result,exp_result,1e-15,"qng(f1) smooth 43pt result") ;
  585.     gsl_test_rel(abserr,exp_abserr,1e-6,"qng(f1) smooth 43pt abserr") ;
  586.     gsl_test_int((int)neval,exp_neval,"qng(f1) smooth 43pt neval") ;  
  587.     gsl_test_int(status,exp_ier,"qng(f1) smooth 43pt status") ;
  588.  
  589.     status = gsl_integration_qng (&f, 1.0, 0.0, 0.0, 1e-9,
  590.                   &result, &abserr, &neval) ;
  591.     gsl_test_rel(result,-exp_result,1e-15,"qng(f1) reverse 43pt result") ;
  592.     gsl_test_rel(abserr,exp_abserr,1e-6,"qng(f1) reverse 43pt abserr") ;
  593.     gsl_test_int((int)neval,exp_neval,"qng(f1) reverse 43pt neval") ;  
  594.     gsl_test_int(status,exp_ier,"qng(f1) reverse 43pt status") ;
  595.   }
  596.  
  597.   {
  598.     int status; size_t neval = 0 ;
  599.     double result = 0, abserr = 0 ;
  600.     double exp_result =-7.238969575482961938E-01;
  601.     double exp_abserr = 1.277676889520056369E-14;
  602.     int exp_neval  =  43;
  603.     int exp_ier    =   0;
  604.  
  605.     double alpha = 1.3 ;
  606.     gsl_function f = make_function(&f3, &alpha);
  607.  
  608.     status = gsl_integration_qng (&f, 0.3, 2.71, 0.0, 1e-12,
  609.                   &result, &abserr, &neval) ;
  610.     gsl_test_rel(result,exp_result,1e-15,"qnq(f3) oscill result") ;
  611.     gsl_test_rel(abserr,exp_abserr,1e-7,"qng(f3) oscill abserr") ;
  612.     gsl_test_int((int)neval,exp_neval,"qng(f3) oscill neval") ;
  613.     gsl_test_int(status,exp_ier,"qng(f3) oscill status") ;
  614.  
  615.     status = gsl_integration_qng (&f, 2.71, 0.3, 0.0, 1e-12,
  616.                   &result, &abserr, &neval) ;
  617.     gsl_test_rel(result,-exp_result,1e-15,"qnq(f3) reverse result") ;
  618.     gsl_test_rel(abserr,exp_abserr,1e-7,"qng(f3) reverse abserr") ;
  619.     gsl_test_int((int)neval,exp_neval,"qng(f3) reverse neval") ;
  620.     gsl_test_int(status,exp_ier,"qng(f3) reverse status") ;
  621.   }
  622.  
  623.   {
  624.     int status = 0; size_t neval = 0 ;
  625.     double result = 0, abserr = 0 ;
  626.  
  627.     double exp_result = 7.716049382716029525E-02;
  628.     double exp_abserr = 8.566535680046930668E-16;
  629.     int exp_neval  =  87;
  630.     int exp_ier    =   0;
  631.  
  632.     double alpha = 2.6 ;
  633.     gsl_function f = make_function(&f1, &alpha);
  634.  
  635.     status = gsl_integration_qng (&f, 0.0, 1.0, 0.0, 1e-13,
  636.                   &result, &abserr, &neval) ;
  637.     gsl_test_rel(result,exp_result,1e-15,"qng(f1) 87pt smooth result") ;
  638.     gsl_test_rel(abserr,exp_abserr,1e-7,"qng(f1) 87pt smooth abserr") ;
  639.     gsl_test_int((int)neval,exp_neval,"qng(f1) 87pt smooth neval") ;  
  640.     gsl_test_int(status,exp_ier,"qng(f1) 87pt smooth status") ;
  641.  
  642.     status = gsl_integration_qng (&f, 1.0, 0.0, 0.0, 1e-13,
  643.                   &result, &abserr, &neval) ;
  644.     gsl_test_rel(result,-exp_result,1e-15,"qng(f1) 87pt reverse result") ;
  645.     gsl_test_rel(abserr,exp_abserr,1e-7,"qng(f1) 87pt reverse abserr") ;
  646.     gsl_test_int((int)neval,exp_neval,"qng(f1) 87pt reverse neval") ;  
  647.     gsl_test_int(status,exp_ier,"qng(f1) 87pt reverse status") ;
  648.   }
  649.  
  650.   {
  651.     int status = 0; size_t neval = 0 ;
  652.     double result = 0, abserr = 0 ;
  653.  
  654.     double exp_result = 3.222948711817264211E+01;
  655.     double exp_abserr = 2.782360287710622870E+01;
  656.     int exp_neval  =  87;
  657.     int exp_ier    =  GSL_ETOL;
  658.  
  659.     double alpha = -0.9 ;
  660.     gsl_function f = make_function(&f1, &alpha);
  661.  
  662.     status = gsl_integration_qng (&f, 0.0, 1.0, 0.0, 1e-3,
  663.                   &result, &abserr, &neval) ;
  664.     gsl_test_rel(result,exp_result,1e-15,"qng(f1) sing beyond 87pt result");
  665.     gsl_test_rel(abserr,exp_abserr,1e-7,"qng(f1) sing beyond 87pt abserr");
  666.     gsl_test_int((int)neval,exp_neval,"qng(f1) sing beyond 87pt neval") ;  
  667.     gsl_test_int(status,exp_ier,"qng(f1) sing beyond 87pt status") ;
  668.  
  669.     status = gsl_integration_qng (&f, 1.0, 0.0, 0.0, 1e-3,
  670.                   &result, &abserr, &neval) ;
  671.     gsl_test_rel(result,-exp_result,1e-15,"qng(f1) reverse beyond 87pt result");
  672.     gsl_test_rel(abserr,exp_abserr,1e-7,"qng(f1) rev beyond 87pt abserr");
  673.     gsl_test_int((int)neval,exp_neval,"qng(f1) rev beyond 87pt neval") ;  
  674.     gsl_test_int(status,exp_ier,"qng(f1) rev beyond 87pt status") ;
  675.   }
  676.  
  677.   /* Test the adaptive integrator QAG */
  678.  
  679.   {
  680.     int status = 0, i; struct counter_params p;
  681.     double result = 0, abserr=0;
  682.  
  683.     gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000) ;
  684.  
  685.     double exp_result = 7.716049382715854665E-02 ;
  686.     double exp_abserr = 6.679384885865053037E-12 ;
  687.     int exp_neval  =     165;
  688.     int exp_ier    =       0;
  689.     int exp_last   =       6;
  690.  
  691.     double a[6] = { 0, 0.5, 0.25, 0.125, 0.0625, 0.03125 } ;
  692.     double b[6] = { 0.03125, 1, 0.5, 0.25, 0.125, 0.0625 } ;
  693.     double r[6] = { 3.966769831709074375E-06, 5.491842501998222409E-02,
  694.             1.909827770934243926E-02, 2.776531175604360531E-03,
  695.             3.280661030752063693E-04, 3.522704932261797744E-05 } ;
  696.     double e[6] = { 6.678528276336181873E-12, 6.097169993333454062E-16,
  697.             2.120334764359736934E-16, 3.082568839745514608E-17,
  698.             3.642265412331439511E-18, 3.910988124757650942E-19 } ;
  699.     int order[6] = { 1, 2, 3, 4, 5, 6 } ;
  700.  
  701.     double alpha = 2.6 ;
  702.     gsl_function f = make_function(&f1, &alpha) ;
  703.  
  704.     gsl_function fc = make_counter(&f, &p) ;
  705.  
  706.     status = gsl_integration_qag (&fc, 0.0, 1.0, 0.0, 1e-10, w->limit,
  707.                   GSL_INTEG_GAUSS15, w,
  708.                   &result, &abserr) ;
  709.  
  710.     gsl_test_rel(result,exp_result,1e-15,"qag(f1) smooth result") ;
  711.     gsl_test_rel(abserr,exp_abserr,1e-6,"qag(f1) smooth abserr") ;
  712.     gsl_test_int((int)(p.neval),exp_neval,"qag(f1) smooth neval") ;  
  713.     gsl_test_int((int)(w->size),exp_last,"qag(f1) smooth last") ;  
  714.     gsl_test_int(status,exp_ier,"qag(f1) smooth status") ;
  715.  
  716.     for (i = 0; i < 6 ; i++) 
  717.     gsl_test_rel(w->alist[i],a[i],1e-15,"qag(f1) smooth alist") ;
  718.  
  719.     for (i = 0; i < 6 ; i++) 
  720.     gsl_test_rel(w->blist[i],b[i],1e-15,"qag(f1) smooth blist") ;
  721.  
  722.     for (i = 0; i < 6 ; i++) 
  723.     gsl_test_rel(w->rlist[i],r[i],1e-15,"qag(f1) smooth rlist") ;
  724.  
  725.     for (i = 0; i < 6 ; i++) 
  726.     gsl_test_rel(w->elist[i],e[i],1e-6,"qag(f1) smooth elist") ;
  727.  
  728.     for (i = 0; i < 6 ; i++) 
  729.     gsl_test_int((int)w->order[i],order[i]-1,"qag(f1) smooth order") ;
  730.  
  731.     p.neval = 0;
  732.  
  733.     status = gsl_integration_qag (&fc, 1.0, 0.0, 0.0, 1e-10, w->limit,
  734.                   GSL_INTEG_GAUSS15, w,
  735.                   &result, &abserr) ;
  736.  
  737.     gsl_test_rel(result,-exp_result,1e-15,"qag(f1) reverse result") ;
  738.     gsl_test_rel(abserr,exp_abserr,1e-6,"qag(f1) reverse abserr") ;
  739.     gsl_test_int((int)(p.neval),exp_neval,"qag(f1) reverse neval") ;  
  740.     gsl_test_int((int)(w->size),exp_last,"qag(f1) reverse last") ;  
  741.     gsl_test_int(status,exp_ier,"qag(f1) reverse status") ;
  742.  
  743.     gsl_integration_workspace_free (w) ;
  744.  
  745.   }
  746.  
  747.   /* Test the same function using an absolute error bound and the
  748.      21-point rule */
  749.  
  750.   {
  751.     int status = 0, i; struct counter_params p;
  752.     double result = 0, abserr=0;
  753.  
  754.     gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000) ;
  755.  
  756.     double exp_result = 7.716049382716050342E-02 ;
  757.     double exp_abserr = 2.227969521869139532E-15 ;
  758.     int exp_neval  =     315;
  759.     int exp_ier    =       0;
  760.     int exp_last   =       8;
  761.  
  762.     double a[8] = { 0, 0.5, 0.25, 0.125, 0.0625, 0.03125, 0.015625,
  763.             0.0078125 } ;
  764.     double b[8] = { 0.0078125, 1, 0.5, 0.25, 0.125, 0.0625, 0.03125,
  765.             0.015625 } ;
  766.     double r[8] = { 3.696942726831556522E-08, 5.491842501998223103E-02,
  767.             1.909827770934243579E-02, 2.776531175604360097E-03,
  768.             3.280661030752062609E-04, 3.522704932261797744E-05,
  769.             3.579060884684503576E-06, 3.507395216921808047E-07 } ;
  770.     double e[8] = { 1.371316364034059572E-15, 6.097169993333454062E-16,
  771.             2.120334764359736441E-16, 3.082568839745514608E-17,
  772.             3.642265412331439511E-18, 3.910988124757650460E-19,
  773.             3.973555800712018091E-20, 3.893990926286736620E-21 } ;
  774.     int order[8] = { 1, 2, 3, 4, 5, 6, 7, 8 } ;
  775.  
  776.     double alpha = 2.6 ;
  777.     gsl_function f = make_function(&f1, &alpha);
  778.  
  779.     gsl_function fc = make_counter(&f, &p) ;
  780.  
  781.     status = gsl_integration_qag (&fc, 0.0, 1.0, 1e-14, 0.0, w->limit,
  782.                   GSL_INTEG_GAUSS21, w,
  783.                   &result, &abserr) ;
  784.  
  785.     gsl_test_rel(result,exp_result,1e-15,"qag(f1,21pt) smooth result") ;
  786.     gsl_test_rel(abserr,exp_abserr,1e-6,"qag(f1,21pt) smooth abserr") ;
  787.     gsl_test_int((int)(p.neval),exp_neval,"qag(f1,21pt) smooth neval") ;  
  788.     gsl_test_int((int)(w->size),exp_last,"qag(f1,21pt) smooth last") ;  
  789.     gsl_test_int(status,exp_ier,"qag(f1,21pt) smooth status") ;
  790.  
  791.     for (i = 0; i < 8 ; i++) 
  792.     gsl_test_rel(w->alist[i],a[i],1e-15,"qag(f1,21pt) smooth alist") ;
  793.  
  794.     for (i = 0; i < 8 ; i++) 
  795.     gsl_test_rel(w->blist[i],b[i],1e-15,"qag(f1,21pt) smooth blist") ;
  796.  
  797.     for (i = 0; i < 8 ; i++) 
  798.     gsl_test_rel(w->rlist[i],r[i],1e-15,"qag(f1,21pt) smooth rlist") ;
  799.  
  800.     for (i = 0; i < 8 ; i++) 
  801.     gsl_test_rel(w->elist[i],e[i],1e-6,"qag(f1,21pt) smooth elist") ;
  802.  
  803.     for (i = 0; i < 8 ; i++) 
  804.     gsl_test_int((int)w->order[i],order[i]-1,"qag(f1,21pt) smooth order");
  805.  
  806.  
  807.     p.neval = 0;
  808.     status = gsl_integration_qag (&fc, 1.0, 0.0, 1e-14, 0.0, w->limit,
  809.                   GSL_INTEG_GAUSS21, w,
  810.                   &result, &abserr) ;
  811.  
  812.     gsl_test_rel(result,-exp_result,1e-15,"qag(f1,21pt) reverse result") ;
  813.     gsl_test_rel(abserr,exp_abserr,1e-6,"qag(f1,21pt) reverse abserr") ;
  814.     gsl_test_int((int)(p.neval),exp_neval,"qag(f1,21pt) reverse neval") ;  
  815.     gsl_test_int((int)(w->size),exp_last,"qag(f1,21pt) reverse last") ;  
  816.     gsl_test_int(status,exp_ier,"qag(f1,21pt) reverse status") ;
  817.  
  818.     gsl_integration_workspace_free (w) ;
  819.  
  820.   }
  821.  
  822.   /* Adaptive integration of an oscillatory function which terminates because
  823.      of roundoff error, uses the 31-pt rule */
  824.  
  825.   {
  826.     int status = 0; struct counter_params p;
  827.     double result = 0, abserr=0;
  828.  
  829.     gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000) ;
  830.  
  831.     double exp_result = -7.238969575482959717E-01;
  832.     double exp_abserr =  1.285805464427459261E-14;
  833.     int exp_neval   =     31;
  834.     int exp_ier     =     GSL_EROUND;
  835.     int exp_last    =     1;
  836.  
  837.     double alpha = 1.3 ;
  838.     gsl_function f = make_function(&f3, &alpha);
  839.  
  840.     gsl_function fc = make_counter(&f, &p) ;
  841.  
  842.     status = gsl_integration_qag (&fc, 0.3, 2.71, 1e-14, 0.0, w->limit, 
  843.                   GSL_INTEG_GAUSS31, w, 
  844.                   &result, &abserr) ;
  845.  
  846.     gsl_test_rel(result,exp_result,1e-15,"qag(f3,31pt) oscill result");
  847.     gsl_test_rel(abserr,exp_abserr,1e-6,"qag(f3,31pt) oscill abserr");
  848.     gsl_test_int((int)(p.neval),exp_neval,"qag(f3,31pt) oscill neval") ;  
  849.     gsl_test_int((int)(w->size),exp_last,"qag(f3,31pt) oscill last") ;  
  850.     gsl_test_int(status,exp_ier,"qag(f3,31pt) oscill status") ;
  851.  
  852.     p.neval = 0;
  853.     status = gsl_integration_qag (&fc, 2.71, 0.3, 1e-14, 0.0, w->limit, 
  854.                   GSL_INTEG_GAUSS31, w, 
  855.                   &result, &abserr) ;
  856.  
  857.     gsl_test_rel(result,-exp_result,1e-15,"qag(f3,31pt) reverse result");
  858.     gsl_test_rel(abserr,exp_abserr,1e-6,"qag(f3,31pt) reverse abserr");
  859.     gsl_test_int((int)(p.neval),exp_neval,"qag(f3,31pt) reverse neval") ;  
  860.     gsl_test_int((int)(w->size),exp_last,"qag(f3,31pt) reverse last") ;  
  861.     gsl_test_int(status,exp_ier,"qag(f3,31pt) reverse status") ;
  862.  
  863.     gsl_integration_workspace_free (w) ;
  864.  
  865.   }
  866.  
  867.   /* Check the singularity detection (singularity at x=-0.1 in this example) */
  868.  
  869.   {
  870.     int status = 0; struct counter_params p;
  871.     double result = 0, abserr=0;
  872.  
  873.     gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000) ;
  874.  
  875.     int exp_neval  =     5151;
  876.     int exp_ier    =     GSL_ESING;
  877.     int exp_last   =     51;
  878.  
  879.     double alpha = 2.0 ;
  880.     gsl_function f = make_function(&f16, &alpha);
  881.  
  882.     gsl_function fc = make_counter(&f, &p) ;
  883.  
  884.     status = gsl_integration_qag (&fc, -1.0, 1.0, 1e-14, 0.0, w->limit,
  885.                   GSL_INTEG_GAUSS51, w, 
  886.                   &result, &abserr) ;
  887.  
  888.     gsl_test_int((int)(p.neval),exp_neval,"qag(f16,51pt) sing neval") ;  
  889.     gsl_test_int((int)(w->size),exp_last,"qag(f16,51pt) sing last") ;  
  890.     gsl_test_int(status,exp_ier,"qag(f16,51pt) sing status") ;
  891.  
  892.     p.neval = 0;
  893.     status = gsl_integration_qag (&fc, 1.0, -1.0, 1e-14, 0.0, w->limit,
  894.                   GSL_INTEG_GAUSS51, w, 
  895.                   &result, &abserr) ;
  896.  
  897.     gsl_test_int((int)(p.neval),exp_neval,"qag(f16,51pt) rev neval") ;  
  898.     gsl_test_int((int)(w->size),exp_last,"qag(f16,51pt) rev last") ;  
  899.     gsl_test_int(status,exp_ier,"qag(f16,51pt) rev status") ;
  900.  
  901.     gsl_integration_workspace_free (w) ;
  902.  
  903.   }
  904.  
  905.   /* Check for hitting the iteration limit */
  906.  
  907.   {
  908.     int status = 0, i; struct counter_params p;
  909.     double result = 0, abserr=0;
  910.  
  911.     gsl_integration_workspace * w = gsl_integration_workspace_alloc (3) ;
  912.  
  913.     double exp_result =  9.565151449233894709 ;
  914.     double exp_abserr =  1.570369823891028460E+01;
  915.     int exp_neval  =     305;
  916.     int exp_ier    =     GSL_EMAXITER;
  917.     int exp_last   =     3;
  918.  
  919.     double a[3] = { -5.000000000000000000E-01,
  920.             0.000000000000000000,
  921.             -1.000000000000000000 } ;
  922.     
  923.     double b[3] = { 0.000000000000000000,
  924.                     1.000000000000000000,
  925.                     -5.000000000000000000E-01 } ;
  926.     
  927.     double r[3] = { 9.460353469435913709,
  928.                     9.090909090909091161E-02,
  929.                     1.388888888888888812E-02 } ;
  930.     
  931.     double e[3] = { 1.570369823891028460E+01,
  932.                     1.009293658750142399E-15,
  933.                     1.541976423090495140E-16 } ;
  934.     
  935.     int order[3] = { 1, 2, 3 } ;
  936.  
  937.     double alpha = 1.0 ;
  938.     gsl_function f = make_function(&f16, &alpha);
  939.     gsl_function fc = make_counter(&f, &p) ;
  940.  
  941.     status = gsl_integration_qag (&fc, -1.0, 1.0, 1e-14, 0.0, w->limit, 
  942.                   GSL_INTEG_GAUSS61, w, 
  943.                   &result, &abserr) ;
  944.  
  945.     gsl_test_rel(result,exp_result,1e-15,"qag(f16,61pt) limit result") ;
  946.     gsl_test_rel(abserr,exp_abserr,1e-6,"qag(f16,61pt) limit abserr") ;
  947.     gsl_test_int((int)(p.neval),exp_neval,"qag(f16,61pt) limit neval") ;  
  948.     gsl_test_int((int)(w->size),exp_last,"qag(f16,61pt) limit last") ;  
  949.     gsl_test_int(status,exp_ier,"qag(f16,61pt) limit status") ;
  950.  
  951.     for (i = 0; i < 3 ; i++) 
  952.     gsl_test_rel(w->alist[i],a[i],1e-15,"qag(f16,61pt) limit alist") ;
  953.  
  954.     for (i = 0; i < 3 ; i++) 
  955.     gsl_test_rel(w->blist[i],b[i],1e-15,"qag(f16,61pt) limit blist") ;
  956.  
  957.     for (i = 0; i < 3 ; i++) 
  958.     gsl_test_rel(w->rlist[i],r[i],1e-15,"qag(f16,61pt) limit rlist") ;
  959.  
  960.     for (i = 0; i < 3 ; i++) 
  961.     gsl_test_rel(w->elist[i],e[i],1e-6,"qag(f16,61pt) limit elist") ;
  962.  
  963.     for (i = 0; i < 3 ; i++) 
  964.     gsl_test_int((int)w->order[i],order[i]-1,"qag(f16,61pt) limit order");
  965.  
  966.     p.neval = 0;
  967.     status = gsl_integration_qag (&fc, 1.0, -1.0, 1e-14, 0.0, w->limit, 
  968.                   GSL_INTEG_GAUSS61, w, 
  969.                   &result, &abserr) ;
  970.  
  971.     gsl_test_rel(result,-exp_result,1e-15,"qag(f16,61pt) reverse result") ;
  972.     gsl_test_rel(abserr,exp_abserr,1e-6,"qag(f16,61pt) reverse abserr") ;
  973.     gsl_test_int((int)(p.neval),exp_neval,"qag(f16,61pt) reverse neval") ;  
  974.     gsl_test_int((int)(w->size),exp_last,"qag(f16,61pt) reverse last") ;  
  975.     gsl_test_int(status,exp_ier,"qag(f16,61pt) reverse status") ;
  976.  
  977.     gsl_integration_workspace_free (w) ;
  978.  
  979.   }
  980.  
  981.   /* Test the adaptive integrator with extrapolation QAGS */
  982.  
  983.   {
  984.     int status = 0, i; struct counter_params p;
  985.     double result = 0, abserr=0;
  986.  
  987.     gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000) ;
  988.  
  989.     double exp_result = 7.716049382715789440E-02 ;
  990.     double exp_abserr = 2.216394961010438404E-12 ;
  991.     int exp_neval  =     189;
  992.     int exp_ier    =       0;
  993.     int exp_last   =       5;
  994.  
  995.     double a[5] = { 0, 0.5, 0.25, 0.125, 0.0625 } ;
  996.     double b[5] = { 0.0625, 1, 0.5, 0.25, 0.125 } ;
  997.     double r[5] = { 3.919381915366914693E-05,
  998.             5.491842501998223103E-02,
  999.             1.909827770934243579E-02,
  1000.             2.776531175604360097E-03,
  1001.             3.280661030752062609E-04 } ;
  1002.     double e[5] = { 2.215538742580964735E-12,
  1003.             6.097169993333454062E-16,
  1004.             2.120334764359736441E-16,
  1005.             3.082568839745514608E-17,
  1006.             3.642265412331439511E-18 } ;
  1007.     int order[5] = { 1, 2, 3, 4, 5 } ;
  1008.  
  1009.     double alpha = 2.6 ;
  1010.     gsl_function f = make_function(&f1, &alpha);
  1011.     gsl_function fc = make_counter(&f, &p) ;
  1012.  
  1013.     status = gsl_integration_qags (&fc, 0.0, 1.0, 0.0, 1e-10, w->limit,
  1014.                    w, 
  1015.                    &result, &abserr) ;
  1016.  
  1017.     gsl_test_rel(result,exp_result,1e-15,"qags(f1) smooth result") ;
  1018.     gsl_test_rel(abserr,exp_abserr,1e-6,"qags(f1) smooth abserr") ;
  1019.     gsl_test_int((int)(p.neval),exp_neval,"qags(f1) smooth neval") ;  
  1020.     gsl_test_int((int)(w->size),exp_last,"qags(f1) smooth last") ;  
  1021.     gsl_test_int(status,exp_ier,"qags(f1) smooth status") ;
  1022.  
  1023.     for (i = 0; i < 5 ; i++) 
  1024.     gsl_test_rel(w->alist[i],a[i],1e-15,"qags(f1) smooth alist") ;
  1025.  
  1026.     for (i = 0; i < 5 ; i++) 
  1027.     gsl_test_rel(w->blist[i],b[i],1e-15,"qags(f1) smooth blist") ;
  1028.  
  1029.     for (i = 0; i < 5 ; i++) 
  1030.     gsl_test_rel(w->rlist[i],r[i],1e-15,"qags(f1) smooth rlist") ;
  1031.  
  1032.     for (i = 0; i < 5 ; i++) 
  1033.     gsl_test_rel(w->elist[i],e[i],1e-6,"qags(f1) smooth elist") ;
  1034.  
  1035.     for (i = 0; i < 5 ; i++) 
  1036.     gsl_test_int((int)w->order[i],order[i]-1,"qags(f1) smooth order") ;
  1037.  
  1038.     p.neval = 0;
  1039.     status = gsl_integration_qags (&fc, 1.0, 0.0, 0.0, 1e-10, w->limit,
  1040.                    w, 
  1041.                    &result, &abserr) ;
  1042.  
  1043.     gsl_test_rel(result,-exp_result,1e-15,"qags(f1) reverse result") ;
  1044.     gsl_test_rel(abserr,exp_abserr,1e-6,"qags(f1) reverse abserr") ;
  1045.     gsl_test_int((int)(p.neval),exp_neval,"qags(f1) reverse neval") ;  
  1046.     gsl_test_int((int)(w->size),exp_last,"qags(f1) reverse last") ;  
  1047.     gsl_test_int(status,exp_ier,"qags(f1) reverse status") ;
  1048.  
  1049.     gsl_integration_workspace_free (w) ;
  1050.  
  1051.   }
  1052.  
  1053.   /* Test f11 using an absolute error bound */
  1054.  
  1055.   {
  1056.     int status = 0, i; struct counter_params p;
  1057.     double result = 0, abserr=0;
  1058.  
  1059.     gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000) ;
  1060.  
  1061.     /* All results are for GSL_IEEE_MODE=double-precision */
  1062.  
  1063.     double exp_result = -5.908755278982136588E+03 ;
  1064.     double exp_abserr = 1.299646281053874554E-10 ;
  1065.     int exp_neval  =     357;
  1066.     int exp_ier    =       0;
  1067.     int exp_last   =       9;
  1068.  
  1069.     double a[9] = { 1.000000000000000000E+00,
  1070.             5.005000000000000000E+02,
  1071.             2.507500000000000000E+02,
  1072.             1.258750000000000000E+02,
  1073.             6.343750000000000000E+01,
  1074.             3.221875000000000000E+01,
  1075.             1.660937500000000000E+01,
  1076.             8.804687500000000000E+00,
  1077.             4.902343750000000000E+00 } ;
  1078.     double b[9] = { 4.902343750000000000E+00,
  1079.             1.000000000000000000E+03,
  1080.             5.005000000000000000E+02,
  1081.             2.507500000000000000E+02,
  1082.             1.258750000000000000E+02,
  1083.             6.343750000000000000E+01,
  1084.             3.221875000000000000E+01,
  1085.             1.660937500000000000E+01,
  1086.             8.804687500000000000E+00 } ;
  1087.     double r[9] = { -3.890977835520834649E+00,
  1088.             -3.297343675805121620E+03,
  1089.             -1.475904154146372775E+03,
  1090.             -6.517404019686431411E+02,
  1091.             -2.829354222635842007E+02,
  1092.             -1.201692001973227519E+02,
  1093.             -4.959999906099650246E+01,
  1094.             -1.971441499411640308E+01,
  1095.             -7.457032710459004399E+00 } ;
  1096.     double e[9] = { 6.448276035006137169E-11,
  1097.             3.660786868980994028E-11,
  1098.             1.638582774073219226E-11,
  1099.             7.235772003440423011E-12,
  1100.             3.141214202790722909E-12,
  1101.             1.334146129098576244E-12,
  1102.             5.506706097890446534E-13,
  1103.             2.188739744348345039E-13,
  1104.             8.278969410534525339E-14 } ;
  1105.     int order[9] = { 1, 2, 3, 4, 5, 6, 7, 8, 9 } ;
  1106.  
  1107.     double alpha = 2.0 ;
  1108.     gsl_function f = make_function(&f11, &alpha);
  1109.     gsl_function fc = make_counter(&f, &p) ;
  1110.  
  1111.     status = gsl_integration_qags (&fc, 1.0, 1000.0, 1e-7, 0.0, w->limit,
  1112.                    w, 
  1113.                    &result, &abserr) ;
  1114.     
  1115.     gsl_test_rel(result,exp_result,1e-15,"qags(f11) smooth result") ;
  1116.     gsl_test_rel(abserr,exp_abserr,1e-3,"qags(f11) smooth abserr") ;
  1117.     gsl_test_int((int)(p.neval),exp_neval,"qags(f11) smooth neval") ;  
  1118.     gsl_test_int((int)(w->size),exp_last,"qags(f11) smooth last") ;  
  1119.     gsl_test_int(status,exp_ier,"qags(f11) smooth status") ;
  1120.  
  1121.     for (i = 0; i < 9 ; i++) 
  1122.     gsl_test_rel(w->alist[i],a[i],1e-15,"qags(f11) smooth alist") ;
  1123.  
  1124.     for (i = 0; i < 9 ; i++) 
  1125.     gsl_test_rel(w->blist[i],b[i],1e-15,"qags(f11) smooth blist") ;
  1126.  
  1127.     for (i = 0; i < 9 ; i++) 
  1128.     gsl_test_rel(w->rlist[i],r[i],1e-15,"qags(f11) smooth rlist") ;
  1129.  
  1130.     for (i = 0; i < 9 ; i++) 
  1131.     gsl_test_rel(w->elist[i],e[i],1e-5,"qags(f11) smooth elist") ;
  1132.  
  1133.     for (i = 0; i < 9 ; i++) 
  1134.     gsl_test_int((int)w->order[i],order[i]-1,"qags(f11) smooth order");
  1135.  
  1136.     p.neval = 0;
  1137.     status = gsl_integration_qags (&fc, 1000.0, 1.0, 1e-7, 0.0, w->limit,
  1138.                    w, 
  1139.                    &result, &abserr) ;
  1140.     
  1141.     gsl_test_rel(result,-exp_result,1e-15,"qags(f11) reverse result") ;
  1142.     gsl_test_rel(abserr,exp_abserr,1e-3,"qags(f11) reverse abserr") ;
  1143.     gsl_test_int((int)(p.neval),exp_neval,"qags(f11) reverse neval") ;  
  1144.     gsl_test_int((int)(w->size),exp_last,"qags(f11) reverse last") ;  
  1145.     gsl_test_int(status,exp_ier,"qags(f11) reverse status") ;
  1146.  
  1147.     gsl_integration_workspace_free (w) ;
  1148.  
  1149.   }
  1150.  
  1151.   /* Test infinite range integral f455 using a relative error bound */
  1152.  
  1153.   {
  1154.     int status = 0, i; struct counter_params p;
  1155.     double result = 0, abserr=0;
  1156.  
  1157.     gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000) ;
  1158.  
  1159.     /* All results are for GSL_IEEE_MODE=double-precision */
  1160.  
  1161.     double exp_result = -3.616892186127022568E-01 ;
  1162.     double exp_abserr = 3.016716913328831851E-06;
  1163.     int exp_neval  =      285;
  1164.     int exp_ier    =        0;
  1165.     int exp_last   =       10;
  1166.  
  1167.     double a[10] = { 9.687500000000000000E-01,
  1168.              0.000000000000000000E+00,
  1169.              5.000000000000000000E-01,
  1170.              2.500000000000000000E-01,
  1171.              7.500000000000000000E-01,
  1172.              1.250000000000000000E-01,
  1173.              8.750000000000000000E-01,
  1174.              6.250000000000000000E-02,
  1175.              9.375000000000000000E-01,
  1176.              3.125000000000000000E-02 } ;
  1177.     double b[10] = { 1.000000000000000000E+00,
  1178.              3.125000000000000000E-02,
  1179.              7.500000000000000000E-01,
  1180.              5.000000000000000000E-01,
  1181.              8.750000000000000000E-01,
  1182.              2.500000000000000000E-01,
  1183.              9.375000000000000000E-01,
  1184.              1.250000000000000000E-01,
  1185.              9.687500000000000000E-01,
  1186.              6.250000000000000000E-02 } ;
  1187.     double r[10] = { -1.390003415539725340E-01,
  1188.              1.429785306003466313E-03,
  1189.              -1.229943369113085765E-02,
  1190.              2.995321156568048898E-03,
  1191.              -4.980050133751051655E-02,
  1192.              2.785385934678596704E-03,
  1193.              -8.653752279614615461E-02,
  1194.              1.736218164975512294E-03,
  1195.              -8.398745675010892142E-02,
  1196.              1.041689192004495576E-03 } ;
  1197.     double e[10] = { 2.395037249893453013E-02,
  1198.              2.161214992172538524E-04,
  1199.              5.720644840858777846E-14,
  1200.              3.325474514168701167E-17,
  1201.              3.147380432198176412E-14,
  1202.              3.092399597147240624E-17,
  1203.              9.607595030230581153E-16,
  1204.              1.927589382528252344E-17,
  1205.              9.324480826368044019E-16,
  1206.              1.156507325466566521E-17 } ;
  1207.     int order[10] = { 1, 2, 3, 5, 7, 9, 4, 6, 8, 10 } ;
  1208.  
  1209.     gsl_function f = make_function(&f455, 0);
  1210.     gsl_function fc = make_counter(&f, &p) ;
  1211.  
  1212.     status = gsl_integration_qagiu (&fc, 0.0, 0.0, 1.0e-3, w->limit,
  1213.                     w, 
  1214.                     &result, &abserr) ;
  1215.     
  1216.     gsl_test_rel(result,exp_result,1e-14,"qagiu(f455) smooth result") ;
  1217.     gsl_test_rel(abserr,exp_abserr,1e-5,"qagiu(f455) smooth abserr") ;
  1218.     gsl_test_int((int)(p.neval),exp_neval,"qagiu(f455) smooth neval") ;  
  1219.     gsl_test_int((int)(w->size),exp_last,"qagiu(f455) smooth last") ;  
  1220.     gsl_test_int(status,exp_ier,"qagiu(f455) smooth status") ;
  1221.  
  1222.     for (i = 0; i < 10 ; i++) 
  1223.     gsl_test_rel(w->alist[i],a[i],1e-15,"qagiu(f455) smooth alist") ;
  1224.  
  1225.     for (i = 0; i < 10 ; i++) 
  1226.     gsl_test_rel(w->blist[i],b[i],1e-15,"qagiu(f455) smooth blist") ;
  1227.  
  1228.     for (i = 0; i < 10 ; i++) 
  1229.     gsl_test_rel(w->rlist[i],r[i],1e-15,"qagiu(f455) smooth rlist") ;
  1230.  
  1231.     for (i = 0; i < 10 ; i++) 
  1232.     gsl_test_rel(w->elist[i],e[i],1e-4,"qagiu(f455) smooth elist") ;
  1233.  
  1234.     for (i = 0; i < 10 ; i++) 
  1235.     gsl_test_int((int)w->order[i],order[i]-1,"qagiu(f455) smooth order");
  1236.  
  1237.     gsl_integration_workspace_free (w) ;
  1238.  
  1239.   }
  1240.  
  1241.   /* Test infinite range integral f15 using a relative error bound */
  1242.  
  1243.   {
  1244.     int status = 0, i; struct counter_params p;
  1245.     double result = 0, abserr=0;
  1246.  
  1247.     gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000) ;
  1248.  
  1249.     /* All results are for GSL_IEEE_MODE=double-precision */
  1250.  
  1251.     double exp_result = 6.553600000000024738E+04;
  1252.     double exp_abserr = 7.121667111456009280E-04;
  1253.     int exp_neval  =      285;
  1254.     int exp_ier    =        0;
  1255.     int exp_last   =       10;
  1256.  
  1257.     double a[10] = { 0.000000000000000000E+00,
  1258.              5.000000000000000000E-01,
  1259.              2.500000000000000000E-01,
  1260.              1.250000000000000000E-01,
  1261.              6.250000000000000000E-02,
  1262.              3.125000000000000000E-02,
  1263.              1.562500000000000000E-02,
  1264.              7.812500000000000000E-03,
  1265.              3.906250000000000000E-03,
  1266.              1.953125000000000000E-03 } ;
  1267.     double b[10] = { 1.953125000000000000E-03,
  1268.              1.000000000000000000E+00,
  1269.              5.000000000000000000E-01,
  1270.              2.500000000000000000E-01,
  1271.              1.250000000000000000E-01,
  1272.              6.250000000000000000E-02,
  1273.              3.125000000000000000E-02,
  1274.              1.562500000000000000E-02,
  1275.              7.812500000000000000E-03,
  1276.              3.906250000000000000E-03 } ;
  1277.     double r[10] = { 1.099297665754340292E+00,
  1278.              3.256176475185617591E-01,
  1279.              8.064694554185326325E+00,
  1280.              8.873128656118993263E+01,
  1281.              6.977679035845269482E+02,
  1282.              4.096981198511257389E+03,
  1283.              1.574317583220441520E+04,
  1284.              2.899418134793237914E+04,
  1285.              1.498314766425578091E+04,
  1286.              9.225251570832365360E+02 } ;
  1287.     double e[10] = { 7.101865971621337814E-04,
  1288.              1.912660677170175771E-08,
  1289.              9.167763417119923333E-08,
  1290.              3.769501719163865578E-07,
  1291.              6.973493131275552509E-07,
  1292.              1.205653952340679711E-07,
  1293.              1.380003928453846583E-07,
  1294.              1.934652413547325474E-07,
  1295.              3.408933028357320364E-07,
  1296.              2.132473175465897029E-09 } ;
  1297.     int order[10] = { 1, 5, 4, 9, 8, 7, 6, 3, 2, 10 } ;
  1298.  
  1299.     double alpha = 5.0;
  1300.  
  1301.     gsl_function f = make_function(&f15, &alpha);
  1302.     gsl_function fc = make_counter(&f, &p) ;
  1303.  
  1304.     status = gsl_integration_qagiu (&fc, 0.0, 0.0, 1.0e-7, w->limit,
  1305.                     w, 
  1306.                     &result, &abserr) ;
  1307.     
  1308.     gsl_test_rel(result,exp_result,1e-14,"qagiu(f15) smooth result") ;
  1309.     gsl_test_rel(abserr,exp_abserr,1e-5,"qagiu(f15) smooth abserr") ;
  1310.     gsl_test_int((int)(p.neval),exp_neval,"qagiu(f15) smooth neval") ;  
  1311.     gsl_test_int((int)(w->size),exp_last,"qagiu(f15) smooth last") ;  
  1312.     gsl_test_int(status,exp_ier,"qagiu(f15) smooth status") ;
  1313.  
  1314.     for (i = 0; i < 10 ; i++) 
  1315.     gsl_test_rel(w->alist[i],a[i],1e-15,"qagiu(f15) smooth alist") ;
  1316.  
  1317.     for (i = 0; i < 10 ; i++) 
  1318.     gsl_test_rel(w->blist[i],b[i],1e-15,"qagiu(f15) smooth blist") ;
  1319.  
  1320.     for (i = 0; i < 10 ; i++) 
  1321.     gsl_test_rel(w->rlist[i],r[i],1e-15,"qagiu(f15) smooth rlist") ;
  1322.  
  1323.     for (i = 0; i < 10 ; i++) 
  1324.     gsl_test_rel(w->elist[i],e[i],1e-4,"qagiu(f15) smooth elist") ;
  1325.  
  1326.     for (i = 0; i < 10 ; i++) 
  1327.     gsl_test_int((int)w->order[i],order[i]-1,"qagiu(f15) smooth order");
  1328.  
  1329.     gsl_integration_workspace_free (w) ;
  1330.  
  1331.   }
  1332.  
  1333.   /* Test infinite range integral f16 using an absolute error bound */
  1334.  
  1335.   {
  1336.     int status = 0, i; struct counter_params p;
  1337.     double result = 0, abserr=0;
  1338.  
  1339.     gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000) ;
  1340.  
  1341.     /* All results are for GSL_IEEE_MODE=double-precision */
  1342.  
  1343.     double exp_result = 1.000000000006713292E-04;
  1344.     double exp_abserr = 3.084062020905636316E-09;
  1345.     int exp_neval  =      165;
  1346.     int exp_ier    =        0;
  1347.     int exp_last   =        6;
  1348.  
  1349.     double a[6] = { 0.000000000000000000E+00,
  1350.             5.000000000000000000E-01,
  1351.             2.500000000000000000E-01,
  1352.             1.250000000000000000E-01,
  1353.             6.250000000000000000E-02,
  1354.             3.125000000000000000E-02 } ;
  1355.     double b[6] = { 3.125000000000000000E-02,
  1356.             1.000000000000000000E+00,
  1357.             5.000000000000000000E-01,
  1358.             2.500000000000000000E-01,
  1359.             1.250000000000000000E-01,
  1360.             6.250000000000000000E-02 } ;
  1361.     double r[6] = { 7.633587786326674618E-05,
  1362.             9.900990099009899620E-07,
  1363.             1.922522349322310737E-06,
  1364.             3.629434715543053753E-06,
  1365.             6.501422186103209199E-06,
  1366.             1.062064387653501389E-05 } ;
  1367.     double e[6] = { 3.084061858351569051E-09,
  1368.             3.112064814755089674E-17,
  1369.             4.543453652226561245E-17,
  1370.             4.908618166361344548E-17,
  1371.             3.014338672269481784E-17,
  1372.             6.795996738013555461E-18 } ;
  1373.     int order[6] = { 1, 4, 3, 2, 5, 6 } ;
  1374.  
  1375.     double alpha = 1.0;
  1376.  
  1377.     gsl_function f = make_function(&f16, &alpha);
  1378.     gsl_function fc = make_counter(&f, &p) ;
  1379.  
  1380.     status = gsl_integration_qagiu (&fc, 99.9, 1.0e-7, 0.0, w->limit,
  1381.                     w, 
  1382.                     &result, &abserr) ;
  1383.     
  1384.     gsl_test_rel(result,exp_result,1e-14,"qagiu(f16) smooth result") ;
  1385.     gsl_test_rel(abserr,exp_abserr,1e-5,"qagiu(f16) smooth abserr") ;
  1386.     gsl_test_int((int)(p.neval),exp_neval,"qagiu(f16) smooth neval") ;  
  1387.     gsl_test_int((int)(w->size),exp_last,"qagiu(f16) smooth last") ;  
  1388.     gsl_test_int(status,exp_ier,"qagiu(f16) smooth status") ;
  1389.  
  1390.     for (i = 0; i < 6 ; i++) 
  1391.     gsl_test_rel(w->alist[i],a[i],1e-15,"qagiu(f16) smooth alist") ;
  1392.  
  1393.     for (i = 0; i < 6 ; i++) 
  1394.     gsl_test_rel(w->blist[i],b[i],1e-15,"qagiu(f16) smooth blist") ;
  1395.  
  1396.     for (i = 0; i < 6 ; i++) 
  1397.     gsl_test_rel(w->rlist[i],r[i],1e-15,"qagiu(f16) smooth rlist") ;
  1398.  
  1399.     for (i = 0; i < 6 ; i++) 
  1400.     gsl_test_rel(w->elist[i],e[i],1e-4,"qagiu(f16) smooth elist") ;
  1401.  
  1402.     for (i = 0; i < 6 ; i++) 
  1403.     gsl_test_int((int)w->order[i],order[i]-1,"qagiu(f16) smooth order");
  1404.  
  1405.     gsl_integration_workspace_free (w) ;
  1406.  
  1407.   }
  1408.  
  1409.   /* Test infinite range integral myfn1 using an absolute error bound */
  1410.  
  1411.   {
  1412.     int status = 0, i; struct counter_params p;
  1413.     double result = 0, abserr=0;
  1414.  
  1415.     gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000) ;
  1416.  
  1417.     /* All results are for GSL_IEEE_MODE=double-precision */
  1418.  
  1419.     double exp_result = 2.275875794468747770E+00;
  1420.     double exp_abserr = 7.436490118267390744E-09;
  1421.     int exp_neval  =      270;
  1422.     int exp_ier    =        0;
  1423.     int exp_last   =        5;
  1424.  
  1425.     double a[5] = { 1.250000000000000000E-01,
  1426.             5.000000000000000000E-01,
  1427.             2.500000000000000000E-01,
  1428.             0.000000000000000000E+00,
  1429.             3.750000000000000000E-01 } ;
  1430.     double b[5] = { 2.500000000000000000E-01,
  1431.             1.000000000000000000E+00,
  1432.             3.750000000000000000E-01,
  1433.             1.250000000000000000E-01,
  1434.             5.000000000000000000E-01 } ;
  1435.     double r[5] = { 4.639317228058405717E-04,
  1436.             1.691664195356748834E+00,
  1437.             1.146307471900291086E-01,
  1438.             4.379392477350953574E-20,
  1439.             4.691169201991640669E-01 } ;
  1440.     double e[5] = { 3.169263960393051137E-09,
  1441.             4.265988974874425043E-09,
  1442.             1.231954072964969637E-12,
  1443.             8.360902986775307673E-20,
  1444.             5.208244060463541433E-15 } ;
  1445.     int order[5] = { 2, 1, 3, 5, 4 } ;
  1446.  
  1447.     gsl_function f = make_function(&myfn1, 0);
  1448.     gsl_function fc = make_counter(&f, &p) ;
  1449.  
  1450.     status = gsl_integration_qagi (&fc, 1.0e-7, 0.0, w->limit,
  1451.                    w, 
  1452.                    &result, &abserr) ;
  1453.     
  1454.     gsl_test_rel(result,exp_result,1e-14,"qagiu(myfn1) smooth result") ;
  1455.     gsl_test_rel(abserr,exp_abserr,1e-5,"qagiu(myfn1) smooth abserr") ;
  1456.     gsl_test_int((int)(p.neval),exp_neval,"qagiu(myfn1) smooth neval") ;  
  1457.     gsl_test_int((int)(w->size),exp_last,"qagiu(myfn1) smooth last") ;  
  1458.     gsl_test_int(status,exp_ier,"qagiu(myfn1) smooth status") ;
  1459.  
  1460.     for (i = 0; i < 5 ; i++) 
  1461.     gsl_test_rel(w->alist[i],a[i],1e-15,"qagiu(myfn1) smooth alist") ;
  1462.  
  1463.     for (i = 0; i < 5 ; i++) 
  1464.     gsl_test_rel(w->blist[i],b[i],1e-15,"qagiu(myfn1) smooth blist") ;
  1465.  
  1466.     for (i = 0; i < 5 ; i++) 
  1467.     gsl_test_rel(w->rlist[i],r[i],1e-14,"qagiu(myfn1) smooth rlist") ;
  1468.  
  1469.     for (i = 0; i < 5 ; i++) 
  1470.     gsl_test_rel(w->elist[i],e[i],1e-4,"qagiu(myfn1) smooth elist") ;
  1471.  
  1472.     for (i = 0; i < 5 ; i++) 
  1473.     gsl_test_int((int)w->order[i],order[i]-1,"qagiu(myfn1) smooth order");
  1474.  
  1475.     gsl_integration_workspace_free (w) ;
  1476.  
  1477.   }
  1478.  
  1479.   /* Test infinite range integral myfn1 using an absolute error bound */
  1480.  
  1481.   {
  1482.     int status = 0, i; struct counter_params p;
  1483.     double result = 0, abserr=0;
  1484.  
  1485.     gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000) ;
  1486.  
  1487.     /* All results are for GSL_IEEE_MODE=double-precision */
  1488.  
  1489.     double exp_result = 2.718281828459044647E+00;
  1490.     double exp_abserr = 1.588185109253204805E-10;
  1491.     int exp_neval  =      135;
  1492.     int exp_ier    =        0;
  1493.     int exp_last   =        5;
  1494.  
  1495.     double a[5] = { 0.000000000000000000E+00,
  1496.             5.000000000000000000E-01,
  1497.             2.500000000000000000E-01,
  1498.             1.250000000000000000E-01,
  1499.             6.250000000000000000E-02 } ;
  1500.     double b[5] = { 6.250000000000000000E-02,
  1501.             1.000000000000000000E+00,
  1502.             5.000000000000000000E-01,
  1503.             2.500000000000000000E-01,
  1504.             1.250000000000000000E-01 } ;
  1505.     double r[5] = { 8.315287189746029816E-07,
  1506.             1.718281828459045091E+00,
  1507.             8.646647167633871867E-01,
  1508.             1.328565310599463256E-01,
  1509.             2.477920647947255521E-03 } ;
  1510.     double e[5] = { 1.533437090413525935E-10,
  1511.             4.117868247943567505E-12,
  1512.             7.802455785301941044E-13,
  1513.             5.395586026138397182E-13,
  1514.             3.713312434866150125E-14 } ;
  1515.     int order[5] = { 1, 2, 3, 4, 5 } ;
  1516.  
  1517.     double alpha = 1.0 ;
  1518.     gsl_function f = make_function(&myfn2, &alpha);
  1519.     gsl_function fc = make_counter(&f, &p) ;
  1520.  
  1521.     status = gsl_integration_qagil (&fc, 1.0, 1.0e-7, 0.0, w->limit,
  1522.                     w, 
  1523.                     &result, &abserr) ;
  1524.     
  1525.     gsl_test_rel(result,exp_result,1e-14,"qagiu(myfn2) smooth result") ;
  1526.     gsl_test_rel(abserr,exp_abserr,1e-5,"qagiu(myfn2) smooth abserr") ;
  1527.     gsl_test_int((int)(p.neval),exp_neval,"qagiu(myfn2) smooth neval") ;  
  1528.     gsl_test_int((int)(w->size),exp_last,"qagiu(myfn2) smooth last") ;  
  1529.     gsl_test_int(status,exp_ier,"qagiu(myfn2) smooth status") ;
  1530.  
  1531.     for (i = 0; i < 5 ; i++) 
  1532.     gsl_test_rel(w->alist[i],a[i],1e-15,"qagiu(myfn2) smooth alist") ;
  1533.  
  1534.     for (i = 0; i < 5 ; i++) 
  1535.     gsl_test_rel(w->blist[i],b[i],1e-15,"qagiu(myfn2) smooth blist") ;
  1536.  
  1537.     for (i = 0; i < 5 ; i++) 
  1538.     gsl_test_rel(w->rlist[i],r[i],1e-14,"qagiu(myfn2) smooth rlist") ;
  1539.  
  1540.     for (i = 0; i < 5 ; i++) 
  1541.     gsl_test_rel(w->elist[i],e[i],1e-4,"qagiu(myfn2) smooth elist") ;
  1542.  
  1543.     for (i = 0; i < 5 ; i++) 
  1544.     gsl_test_int((int)w->order[i],order[i]-1,"qagiu(myfn2) smooth order");
  1545.  
  1546.     gsl_integration_workspace_free (w) ;
  1547.  
  1548.   }
  1549.  
  1550.   /* Test integral f454 with integrable singular points */
  1551.  
  1552.   {
  1553.     int status = 0, i; struct counter_params p;
  1554.     double result = 0, abserr=0;
  1555.  
  1556.     gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000) ;
  1557.  
  1558.     /* All results are for GSL_IEEE_MODE=double-precision */
  1559.  
  1560.     double exp_result = 5.274080611672716401E+01;
  1561.     double exp_abserr = 1.755703848687062418E-04;
  1562.     int exp_neval  =        777;
  1563.     int exp_ier    =          0;
  1564.     int exp_last   =         20;
  1565.  
  1566.     double a[20] = { 9.687500000000000000E-01,
  1567.             1.401269388548935790E+00,
  1568.             1.414213562373095145E+00,
  1569.             1.000000000000000000E+00,
  1570.             0.000000000000000000E+00,
  1571.             2.207106781186547462E+00,
  1572.             1.810660171779821415E+00,
  1573.             1.207106781186547462E+00,
  1574.             5.000000000000000000E-01,
  1575.             1.103553390593273731E+00,
  1576.             1.612436867076458391E+00,
  1577.             1.310660171779821415E+00,
  1578.             7.500000000000000000E-01,
  1579.             1.051776695296636976E+00,
  1580.             1.513325214724776657E+00,
  1581.             1.362436867076458391E+00,
  1582.             8.750000000000000000E-01,
  1583.             1.463769388548935790E+00,
  1584.             1.388325214724776657E+00,
  1585.             9.375000000000000000E-01} ;
  1586.     double b[20] = { 1.000000000000000000E+00,
  1587.              1.414213562373095145E+00,
  1588.              1.463769388548935790E+00,
  1589.              1.051776695296636976E+00,
  1590.              5.000000000000000000E-01,
  1591.              3.000000000000000000E+00,
  1592.              2.207106781186547462E+00,
  1593.              1.310660171779821415E+00,
  1594.              7.500000000000000000E-01,
  1595.              1.207106781186547462E+00,
  1596.              1.810660171779821415E+00,
  1597.              1.362436867076458391E+00,
  1598.              8.750000000000000000E-01,
  1599.              1.103553390593273731E+00,
  1600.              1.612436867076458391E+00,
  1601.              1.388325214724776657E+00,
  1602.              9.375000000000000000E-01,
  1603.              1.513325214724776657E+00,
  1604.              1.401269388548935790E+00,
  1605.              9.687500000000000000E-01} ;
  1606.     double r[20] = { -1.125078814079027711E-01,
  1607.              -1.565132123531515207E-01,
  1608.              -4.225328513207429193E-01,
  1609.              -1.830392049835374568E-01,
  1610.              6.575875041899758092E-03,
  1611.              4.873920540843067783E+01,
  1612.              6.032891565603589079E+00,
  1613.              -2.991531901645863023E-01,
  1614.              -7.326282608704996063E-03,
  1615.              -2.431894410706912923E-01,
  1616.              5.911661670635662835E-01,
  1617.              -2.236786562536174916E-01,
  1618.              -5.647871991778510847E-02,
  1619.              -1.305470403178642658E-01,
  1620.              -1.721363984401322045E-01,
  1621.              -1.589345454585119055E-01,
  1622.              -7.406626263352669715E-02,
  1623.              -2.208730668000830344E-01,
  1624.              -1.048692749517999567E-01,
  1625.              -6.302287584527696551E-02} ;
  1626.     double e[20] = { 2.506431410088378817E-02,
  1627.              2.730454695485963826E-02,
  1628.              1.017446081816190118E-01,
  1629.              3.252808038935910834E-02,
  1630.              7.300687878575027348E-17,
  1631.              5.411138804637469780E-13,
  1632.              6.697855121200013106E-14,
  1633.              3.321267596107916554E-15,
  1634.              1.417509685426979386E-16,
  1635.              2.699945168224041491E-15,
  1636.              6.573952690524728748E-15,
  1637.              2.483331942899818875E-15,
  1638.              6.270397525408045936E-16,
  1639.              1.449363299575615261E-15,
  1640.              1.911097929242846383E-15,
  1641.              1.764527917763735212E-15,
  1642.              8.223007012367522077E-16,
  1643.              2.452183642810224359E-15,
  1644.              1.164282836272345215E-15,
  1645.              6.996944784151910810E-16} ;
  1646.     int order[20] = { 3, 4, 2, 1, 6, 7, 11, 8, 10, 12, 18,
  1647.              15, 16, 14, 19, 17, 20, 13, 9, 5 } ;
  1648.  
  1649.     gsl_function f = make_function(&f454, 0);
  1650.     gsl_function fc = make_counter(&f, &p) ;
  1651.  
  1652.     double pts[4] ;
  1653.  
  1654.     pts[0] = 0.0;
  1655.     pts[1] = 1.0;
  1656.     pts[2] = sqrt(2.0);
  1657.     pts[3] = 3.0;
  1658.  
  1659.     status = gsl_integration_qagp (&fc, pts, 4,
  1660.                    0.0, 1.0e-3, w->limit,
  1661.                    w, 
  1662.                    &result, &abserr) ;
  1663.     
  1664.     gsl_test_rel(result,exp_result,1e-14,"qagp(f454) singular result") ;
  1665.     gsl_test_rel(abserr,exp_abserr,1e-5,"qagp(f454) singular abserr") ;
  1666.     gsl_test_int((int)(p.neval),exp_neval,"qagp(f454) singular neval") ;  
  1667.     gsl_test_int((int)(w->size),exp_last,"qagp(f454) singular last") ;  
  1668.     gsl_test_int(status,exp_ier,"qagp(f454) singular status") ;
  1669.  
  1670.     for (i = 0; i < 20 ; i++) 
  1671.     gsl_test_rel(w->alist[i],a[i],1e-15,"qagp(f454) singular alist") ;
  1672.  
  1673.     for (i = 0; i < 20 ; i++) 
  1674.     gsl_test_rel(w->blist[i],b[i],1e-15,"qagp(f454) singular blist") ;
  1675.  
  1676.     for (i = 0; i < 20 ; i++) 
  1677.     gsl_test_rel(w->rlist[i],r[i],1e-14,"qagp(f454) singular rlist") ;
  1678.  
  1679.     for (i = 0; i < 20 ; i++) 
  1680.     gsl_test_rel(w->elist[i],e[i],1e-4,"qagp(f454) singular elist") ;
  1681.  
  1682.     for (i = 0; i < 20 ; i++) 
  1683.     gsl_test_int((int)w->order[i],order[i]-1,"qagp(f454) singular order");
  1684.  
  1685.     gsl_integration_workspace_free (w) ;
  1686.  
  1687.   }
  1688.  
  1689.  
  1690.   /* Test cauchy integration using a relative error bound */
  1691.  
  1692.   {
  1693.     int status = 0, i; struct counter_params p;
  1694.     double result = 0, abserr=0;
  1695.  
  1696.     gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000) ;
  1697.  
  1698.     /* All results are for GSL_IEEE_MODE=double-precision */
  1699.  
  1700.     double exp_result = -8.994400695837000137E-02;
  1701.     double exp_abserr =  1.185290176227023727E-06;
  1702.     int exp_neval  =      215;
  1703.     int exp_ier    =        0;
  1704.     int exp_last   =        6;
  1705.  
  1706.     double a[6] = { -1.000000000000000000E+00,
  1707.             2.500000000000000000E+00,
  1708.             1.250000000000000000E+00,
  1709.             6.250000000000000000E-01,
  1710.             -5.000000000000000000E-01,
  1711.             -7.500000000000000000E-01} ;
  1712.     double b[6] = { -7.500000000000000000E-01,
  1713.             5.000000000000000000E+00,
  1714.             2.500000000000000000E+00,
  1715.             1.250000000000000000E+00,
  1716.             6.250000000000000000E-01,
  1717.             -5.000000000000000000E-01} ;
  1718.     double r[6] = { -1.234231128040012976E-01,
  1719.             3.579970394639702888E-03,
  1720.             2.249831615049339983E-02,
  1721.             7.214232992127905808E-02,
  1722.             2.079093855884046535E-02,
  1723.             -8.553244917962132821E-02} ;
  1724.     double e[6] = { 1.172832717970022565E-06,
  1725.             9.018232896137375412E-13,
  1726.             1.815172652101790755E-12,
  1727.             1.006998195150956048E-13,
  1728.             1.245463873006391609E-08,
  1729.             1.833082948207153514E-15 } ;
  1730.     int order[6] = { 1, 5, 3, 2, 4, 6 } ;
  1731.  
  1732.     double alpha = 1.0 ;
  1733.     gsl_function f = make_function(&f459, &alpha);
  1734.     gsl_function fc = make_counter(&f, &p) ;
  1735.  
  1736.     status = gsl_integration_qawc (&fc, -1.0, 5.0, 0.0, 0.0, 1.0e-3, w->limit,
  1737.                    w, 
  1738.                    &result, &abserr) ;
  1739.     
  1740.     gsl_test_rel(result,exp_result,1e-14,"qawc(f459) result") ;
  1741.     gsl_test_rel(abserr,exp_abserr,1e-6,"qawc(f459) abserr") ;
  1742.     gsl_test_int((int)(p.neval),exp_neval,"qawc(f459) neval") ;  
  1743.     gsl_test_int((int)(w->size),exp_last,"qawc(f459) last") ;  
  1744.     gsl_test_int(status,exp_ier,"qawc(f459) status") ;
  1745.  
  1746.     for (i = 0; i < 6 ; i++) 
  1747.     gsl_test_rel(w->alist[i],a[i],1e-15,"qawc(f459) alist") ;
  1748.  
  1749.     for (i = 0; i < 6 ; i++) 
  1750.     gsl_test_rel(w->blist[i],b[i],1e-15,"qawc(f459) blist") ;
  1751.  
  1752.     for (i = 0; i < 6 ; i++) 
  1753.     gsl_test_rel(w->rlist[i],r[i],1e-14,"qawc(f459) rlist") ;
  1754.  
  1755.     for (i = 0; i < 6 ; i++) 
  1756.     gsl_test_rel(w->elist[i],e[i],1e-5,"qawc(f459) elist") ;
  1757.  
  1758.     for (i = 0; i < 6 ; i++) 
  1759.     gsl_test_int((int)w->order[i],order[i]-1,"qawc(f459) order");
  1760.  
  1761.     p.neval = 0;
  1762.     status = gsl_integration_qawc (&fc, 5.0, -1.0, 0.0, 0.0, 1.0e-3, w->limit,
  1763.                    w, 
  1764.                    &result, &abserr) ;
  1765.     
  1766.     gsl_test_rel(result,-exp_result,1e-14,"qawc(f459) rev result") ;
  1767.     gsl_test_rel(abserr,exp_abserr,1e-6,"qawc(f459) rev abserr") ;
  1768.     gsl_test_int((int)(p.neval),exp_neval,"qawc(f459) rev neval") ;  
  1769.     gsl_test_int((int)(w->size),exp_last,"qawc(f459) rev last") ;  
  1770.     gsl_test_int(status,exp_ier,"qawc(f459) rev status") ;
  1771.  
  1772.     gsl_integration_workspace_free (w) ;
  1773.  
  1774.   }
  1775.  
  1776.   /* Test QAWS singular integration using a relative error bound */
  1777.  
  1778.   {
  1779.     int status = 0, i; struct counter_params p;
  1780.     double result = 0, abserr=0;
  1781.  
  1782.     gsl_integration_qaws_table * t 
  1783.       = gsl_integration_qaws_table_alloc (0.0, 0.0, 1, 0);
  1784.  
  1785.     gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000) ;
  1786.  
  1787.     /* All results are for GSL_IEEE_MODE=double-precision */
  1788.  
  1789.     double exp_result = -1.892751853489401670E-01;
  1790.     double exp_abserr = 1.129133712015747658E-08;
  1791.     int exp_neval  =      280;
  1792.     int exp_ier    =        0;
  1793.     int exp_last   =        8;
  1794.  
  1795.     double a[8] = { 0.000000000000000000E+00,
  1796.             5.000000000000000000E-01,
  1797.             2.500000000000000000E-01,
  1798.             1.250000000000000000E-01,
  1799.             6.250000000000000000E-02,
  1800.             3.125000000000000000E-02,
  1801.             1.562500000000000000E-02,
  1802.             7.812500000000000000E-03} ;
  1803.     double b[8] = { 7.812500000000000000E-03,
  1804.             1.000000000000000000E+00,
  1805.             5.000000000000000000E-01,
  1806.             2.500000000000000000E-01,
  1807.             1.250000000000000000E-01,
  1808.             6.250000000000000000E-02,
  1809.             3.125000000000000000E-02,
  1810.             1.562500000000000000E-02} ;
  1811.     double r[8] = { -4.126317299834445824E-05,
  1812.             -1.076283950172247789E-01,
  1813.             -6.240573216173390947E-02,
  1814.             -1.456169844189576269E-02,
  1815.             -3.408925115926728436E-03,
  1816.             -8.914083918175634211E-04,
  1817.             -2.574191402137795482E-04,
  1818.             -8.034390712936630608E-05} ;
  1819.     double e[8] = { 1.129099387465713953E-08,
  1820.             3.423394967694403596E-13,
  1821.             6.928428071454762659E-16,
  1822.             1.616673288784094320E-16,
  1823.             3.784667152924835070E-17,
  1824.             9.896621209399419425E-18,
  1825.             2.857926564445496100E-18,
  1826.             8.919965558336773736E-19} ;
  1827.     int order[8] = { 1, 2, 3, 4, 5, 6, 7, 8 } ;
  1828.  
  1829.     double alpha = 1.0 ;
  1830.     gsl_function f = make_function(&f458, &alpha);
  1831.     gsl_function fc = make_counter(&f, &p) ;
  1832.  
  1833.     status = gsl_integration_qaws (&fc, 0.0, 1.0, t, 0.0, 1.0e-7, w->limit,
  1834.                    w, 
  1835.                    &result, &abserr) ;
  1836.     
  1837.     gsl_test_rel(result,exp_result,1e-14,"qaws(f458) ln(x-a) result") ;
  1838.     gsl_test_rel(abserr,exp_abserr,1e-6,"qaws(f458) ln(x-a) abserr") ;
  1839.     gsl_test_int((int)(p.neval),exp_neval,"qaws(f458) ln(x-a) neval") ;  
  1840.     gsl_test_int((int)(w->size),exp_last,"qaws(f458) ln(x-a) last") ;  
  1841.     gsl_test_int(status,exp_ier,"qaws(f458) ln(x-a) status") ;
  1842.  
  1843.     for (i = 0; i < 6 ; i++) 
  1844.     gsl_test_rel(w->alist[i],a[i],1e-15,"qaws(f458) ln(x-a) alist") ;
  1845.  
  1846.     for (i = 0; i < 6 ; i++) 
  1847.     gsl_test_rel(w->blist[i],b[i],1e-15,"qaws(f458) ln(x-a) blist") ;
  1848.  
  1849.     for (i = 0; i < 6 ; i++) 
  1850.     gsl_test_rel(w->rlist[i],r[i],1e-14,"qaws(f458) ln(x-a) rlist") ;
  1851.  
  1852.     for (i = 0; i < 6 ; i++) 
  1853.     gsl_test_rel(w->elist[i],e[i],1e-4,"qaws(f458) ln(x-a) elist") ;
  1854.  
  1855.     for (i = 0; i < 6 ; i++) 
  1856.     gsl_test_int((int)w->order[i],order[i]-1,"qaws(f458) ln(x-a) order");
  1857.     
  1858.     /* Test without logs */
  1859.     
  1860.     gsl_integration_qaws_table_set (t, -0.5, -0.3, 0, 0);
  1861.     
  1862.     status = gsl_integration_qaws (&fc, 0.0, 1.0, t, 0.0, 1.0e-7, w->limit,
  1863.                    w, &result, &abserr) ;
  1864.  
  1865.     exp_result = 9.896686656601706433E-01;
  1866.     exp_abserr = 5.888032513201251628E-08;
  1867.  
  1868.     gsl_test_rel(result,exp_result,1e-14,"qaws(f458) AB result") ;
  1869.     gsl_test_rel(abserr,exp_abserr,1e-6,"qaws(f458) AB abserr") ;
  1870.  
  1871.     /* Test with ln(x - a) */
  1872.  
  1873.     gsl_integration_qaws_table_set (t, -0.5, -0.3, 1, 0);
  1874.     
  1875.     status = gsl_integration_qaws (&fc, 0.0, 1.0, t, 0.0, 1.0e-7, w->limit,
  1876.                    w, &result, &abserr) ;
  1877.  
  1878.     exp_result = -3.636679470586539620E-01;
  1879.     exp_abserr = 2.851348775257054093E-08;
  1880.  
  1881.     gsl_test_rel(result,exp_result,1e-14,"qaws(f458) AB ln(x-a) result") ;
  1882.     gsl_test_rel(abserr,exp_abserr,1e-6,"qaws(f458) AB ln(x-a) abserr") ;
  1883.  
  1884.     /* Test with ln(b - x) */
  1885.  
  1886.     gsl_integration_qaws_table_set (t, -0.5, -0.3, 0, 1);
  1887.     
  1888.     status = gsl_integration_qaws (&fc, 0.0, 1.0, t, 0.0, 1.0e-7, w->limit,
  1889.                    w, &result, &abserr) ;
  1890.  
  1891.     exp_result = -1.911489253363409802E+00;
  1892.     exp_abserr = 9.854016753016499034E-09;
  1893.  
  1894.     gsl_test_rel(result,exp_result,1e-14,"qaws(f458) AB ln(b-x) result") ;
  1895.     gsl_test_rel(abserr,exp_abserr,1e-6,"qaws(f458) AB ln(b-x) abserr") ;
  1896.  
  1897.     /* Test with ln(x - a) ln(b - x) */
  1898.  
  1899.     gsl_integration_qaws_table_set (t, -0.5, -0.3, 1, 1);
  1900.     
  1901.     status = gsl_integration_qaws (&fc, 0.0, 1.0, t, 0.0, 1.0e-7, w->limit,
  1902.                    w, &result, &abserr) ;
  1903.  
  1904.     exp_result = 3.159922862811048172E-01;
  1905.     exp_abserr = 2.336183482198144595E-08;
  1906.  
  1907.     gsl_test_rel(result,exp_result,1e-14,"qaws(f458) AB ln(x-a)ln(b-x) result") ;
  1908.     gsl_test_rel(abserr,exp_abserr,1e-6,"qaws(f458) AB ln(x-a)ln(b-x) abserr") ;
  1909.  
  1910.     gsl_integration_workspace_free (w) ;
  1911.     gsl_integration_qaws_table_free (t) ;
  1912.  
  1913.   }
  1914.  
  1915.  
  1916.   /* Test oscillatory integration using a relative error bound */
  1917.  
  1918.   {
  1919.     int status = 0, i; struct counter_params p;
  1920.     double result = 0, abserr=0;
  1921.  
  1922.     gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000) ;
  1923.     gsl_integration_qawo_table * wo 
  1924.       = gsl_integration_qawo_table_alloc (10.0 * M_PI, 1.0,
  1925.                           GSL_INTEG_SINE, 1000) ;
  1926.  
  1927.     /* All results are for GSL_IEEE_MODE=double-precision */
  1928.  
  1929.     double exp_result = -1.281368483991674190E-01;
  1930.     double exp_abserr =  6.875028324415666248E-12;
  1931.     int exp_neval  =      305;
  1932.     int exp_ier    =        0;
  1933.     int exp_last   =        9;
  1934.  
  1935.     double a[9] = { 0.000000000000000000E+00,
  1936.             5.000000000000000000E-01,
  1937.             2.500000000000000000E-01,
  1938.             1.250000000000000000E-01,
  1939.             6.250000000000000000E-02,
  1940.             3.125000000000000000E-02,
  1941.             1.562500000000000000E-02,
  1942.             7.812500000000000000E-03,
  1943.             3.906250000000000000E-03 } ;
  1944.     double b[9] = { 3.906250000000000000E-03,
  1945.             1.000000000000000000E+00,
  1946.             5.000000000000000000E-01,
  1947.             2.500000000000000000E-01,
  1948.             1.250000000000000000E-01,
  1949.             6.250000000000000000E-02,
  1950.             3.125000000000000000E-02,
  1951.             1.562500000000000000E-02,
  1952.             7.812500000000000000E-03 } ;
  1953.     double r[9] = { -1.447193692377651136E-03,
  1954.             2.190541162282139478E-02,
  1955.             -2.587726479625663753E-02,
  1956.             5.483209176363500886E-02,
  1957.             -3.081695575172510582E-02,
  1958.             -9.178321994387816929E-02,
  1959.             -3.886716016498160953E-02,
  1960.             -1.242306301902117854E-02,
  1961.             -3.659495117871544145E-03} ;
  1962.     double e[9] = { 8.326506625798146465E-07,
  1963.             1.302638552580516100E-13,
  1964.             7.259224351945759794E-15,
  1965.             1.249770395036711102E-14,
  1966.             7.832180081562836579E-16,
  1967.             1.018998440559284116E-15,
  1968.             4.315121611695628020E-16,
  1969.             1.379237060008662177E-16,
  1970.             4.062855738364339357E-17 } ;
  1971.     int order[9] = { 1, 2, 4, 3, 6, 5, 7, 8, 9 } ;
  1972.  
  1973.     double alpha = 1.0 ;
  1974.     gsl_function f = make_function(&f456, &alpha);
  1975.     gsl_function fc = make_counter(&f, &p) ;
  1976.  
  1977.     status = gsl_integration_qawo (&fc, 0.0, 0.0, 1e-7, w->limit,
  1978.                    w, wo, &result, &abserr) ;
  1979.     
  1980.     gsl_test_rel(result,exp_result,1e-14,"qawo(f456) result") ;
  1981.     gsl_test_rel(abserr,exp_abserr,1e-3,"qawo(f456) abserr") ;
  1982.     gsl_test_int((int)(p.neval),exp_neval,"qawo(f456) neval") ;  
  1983.     gsl_test_int((int)(w->size),exp_last,"qawo(f456) last") ;  
  1984.     gsl_test_int(status,exp_ier,"qawo(f456) status") ;
  1985.  
  1986.     for (i = 0; i < 9 ; i++) 
  1987.     gsl_test_rel(w->alist[i],a[i],1e-15,"qawo(f456) alist") ;
  1988.  
  1989.     for (i = 0; i < 9 ; i++) 
  1990.     gsl_test_rel(w->blist[i],b[i],1e-15,"qawo(f456) blist") ;
  1991.  
  1992.     for (i = 0; i < 9 ; i++) 
  1993.     gsl_test_rel(w->rlist[i],r[i],1e-14,"qawo(f456) rlist") ;
  1994.  
  1995.     for (i = 0; i < 9 ; i++) 
  1996.     gsl_test_rel(w->elist[i],e[i],1e-3,"qawo(f456) elist") ;
  1997.  
  1998.     for (i = 0; i < 9 ; i++) 
  1999.     gsl_test_int((int)w->order[i],order[i]-1,"qawo(f456) order");
  2000.  
  2001.  
  2002.     /* In reverse, flip limit and sign of length */
  2003.  
  2004.     gsl_integration_qawo_table_set_length (wo, -1.0);
  2005.  
  2006.     p.neval = 0; 
  2007.     status = gsl_integration_qawo (&fc, 1.0, 0.0, 1e-7, w->limit,
  2008.                    w, wo, &result, &abserr) ;
  2009.     
  2010.     gsl_test_rel(result,-exp_result,1e-14,"qawo(f456) rev result") ;
  2011.     gsl_test_rel(abserr,exp_abserr,1e-3,"qawo(f456) rev abserr") ;
  2012.     gsl_test_int((int)(p.neval),exp_neval,"qawo(f456) rev neval") ;  
  2013.     gsl_test_int((int)(w->size),exp_last,"qawo(f456) rev last") ;  
  2014.     gsl_test_int(status,exp_ier,"qawo(f456) rev status") ;
  2015.  
  2016.  
  2017.     gsl_integration_qawo_table_free (wo) ;
  2018.     gsl_integration_workspace_free (w) ;
  2019.  
  2020.   }
  2021.  
  2022.   /* Test fourier integration using an absolute error bound */
  2023.  
  2024.   {
  2025.     int status = 0, i; struct counter_params p;
  2026.     double result = 0, abserr=0;
  2027.  
  2028.     gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000) ;
  2029.     gsl_integration_workspace * wc = gsl_integration_workspace_alloc (1000) ;
  2030.     gsl_integration_qawo_table * wo 
  2031.       = gsl_integration_qawo_table_alloc (M_PI / 2.0, 1.0,
  2032.                           GSL_INTEG_COSINE, 1000) ;
  2033.  
  2034.     /* All results are for GSL_IEEE_MODE=double-precision */
  2035.  
  2036.     double exp_result = 9.999999999279802765E-01;
  2037.     double exp_abserr = 1.556289974669056164E-08;
  2038.     int exp_neval  =      590;
  2039.     int exp_ier    =        0;
  2040.     int exp_last   =       12;
  2041.  
  2042.     double r[12] = { 1.013283128125232802E+00,
  2043.             -1.810857954748607349E-02,
  2044.             7.466754034900931897E-03,
  2045.             -4.360312526786496237E-03,
  2046.             2.950184068216192904E-03,
  2047.             -2.168238443073697373E-03,
  2048.             1.680910783140869081E-03,
  2049.             -1.352797860944863345E-03,
  2050.             1.119354921991485901E-03,
  2051.             -9.462367583691360827E-04,
  2052.             8.136341270731781887E-04,
  2053.             -7.093931338504278145E-04 } ;
  2054.     double e[12] = { 1.224798040766472695E-12,
  2055.             1.396565155187268456E-13,
  2056.             1.053844511655910310E-16,
  2057.             6.505213034913026604E-19,
  2058.             7.155734338404329264E-18,
  2059.             1.105886215935214523E-17,
  2060.             9.757819552369539906E-18,
  2061.             5.854691731421723944E-18,
  2062.             4.553649124439220312E-18,
  2063.             7.643625316022806260E-18,
  2064.             2.439454888092388058E-17,
  2065.             2.130457268934021451E-17 } ;
  2066.  
  2067.     double alpha = 1.0 ;
  2068.     gsl_function f = make_function(&f457, &alpha);
  2069.     gsl_function fc = make_counter(&f, &p) ;
  2070.  
  2071.     status = gsl_integration_qawf (&fc, 0.0, 1e-7, w->limit,
  2072.                    w, wc, wo, &result, &abserr) ;
  2073.     
  2074.     gsl_test_rel(result,exp_result,1e-14,"qawf(f457) result") ;
  2075.     gsl_test_rel(abserr,exp_abserr,1e-3,"qawf(f457) abserr") ;
  2076.     gsl_test_int((int)(p.neval),exp_neval,"qawf(f457) neval") ;  
  2077.     gsl_test_int((int)(w->size),exp_last,"qawf(f457) last") ;  
  2078.     gsl_test_int(status,exp_ier,"qawf(f457) status") ;
  2079.  
  2080.     for (i = 0; i < 9 ; i++) 
  2081.     gsl_test_rel(w->rlist[i],r[i],1e-12,"qawf(f457) rlist") ;
  2082.  
  2083.     /* We can only get within two orders of magnitude on the error
  2084.        here, which is very sensitive to the floating point precision */
  2085.  
  2086.     for (i = 0; i < 9 ; i++) 
  2087.     gsl_test_rel(w->elist[i],e[i],50.0,"qawf(f457) elist") ;
  2088.  
  2089.  
  2090.     gsl_integration_qawo_table_free (wo) ;
  2091.     gsl_integration_workspace_free (wc) ;
  2092.     gsl_integration_workspace_free (w) ;
  2093.  
  2094.   }
  2095.  
  2096.   exit (gsl_test_summary());
  2097. } 
  2098.  
  2099. void
  2100. my_error_handler (const char *reason, const char *file, int line, int err)
  2101. {
  2102.   if (0) printf ("(caught [%s:%d: %s (%d)])\n", file, line, reason, err) ;
  2103. }
  2104.  
  2105.  
  2106.